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NOMENCLATURE 


C  =  physical  damping  matrix 
M  =  physical  mass  matrix 
K  =  physical  stiffness  matrix 
C  =  modal  damping  matrix 
M  =  modal  mass  matrix,  identity  matrix  [I] 

K  =  modal  stiffness  matrix 
T  =  kinetic  energy 
D  =damping  dissipation  function 
U  =  potential  energy 
X  =  lagrange  multiplier 
®  =  mass  normalized  modal  matrix 
(f>  =  mode  shape 

P  =  another  way  to  represent  mode  shape,  individual  vector,  or  element  of  mode 
shape 

co  =  natural  frequency  Hz 

P  =  modal  damping 

Qd  =  undamped  natural  frequency 

a  =  Cm  ,  modal  damping  times  natural  frequency 

x  =  physical  coordinate 

q  =  modal  coordinate 

[  I  ]  =  identity  matrix 

ix 


/  =  interconnecting  condition  between  substructure 


t  =  time 


Subscripts 
j  =  mode  index 

r  =  number  of  interconnecting  region 
i  =  element  of  substructure 

Superscript 

n  =  substructure  designator 
’  =  transpose 
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I.  INTRODUCTION 


Structural  synthesis  is  the  process  of  calculating  the  response  of  an  assemblage  of 
two  or  more  substructures  where  each  substructure  is  represented  by  impulse  response 
functions  calculated  at  the  connection  degrees-of-freedom  only.  In  addition,  since  it  is 
often  difficult  or  even  impossible  to  conduct  vibrations  tests  on  the  system,  it  is  often 
desired  to  extract  system’s  dynamic  response  by  breaking  the  system  into  substructures 
and  conduct  vibrations  test  on  these  individual  substructures.  The  dynamic  responses 
from  the  substructures  are  then  synthesized  to  provide  the  dynamic  response  for  the 
system.  In  this  thesis  we  will  exploit  the  concept  of  substructure  coupling  which  makes 
use  of  linear  homogenous  equations  of  constraints;  such  a  method  enables  us  to  generate 
the  system  dynamic  response  by  making  use  of  substructure  modal  properties 
(transformed  physical  properties  to  modal  properties)  and  coupling  constraints.  The 
method  that  will  be  implemented  to  calculate  the  system  response  is  based  on  the  zero- 
eigenvalues  theorem  developed  by  Walton  and  Steeves  [16]. 

With  the  ability  to  calculate  system  dynamic  response  we  are  now  able  to  go 
further  and  investigate  a  method  of  determining  bounds  for  substructure  damping  that 
would  satisfy  prescribed  bounds  for  system  damping  given  only  substructure  physical 
properties  such  as  mass  and  stiffness  and  coupling  constraints.  In  support  of  the  task  in 
hand,  we  will  also  make  use  of  other  concepts  such  as  Lagrange’s  equation  to  attain 
equations  of  motion,  Lagrange  multipliers  and  holonomic  constraints  to  provide 
connectivity  for  each  substructure,  orthogonal  complement  and  zero-eigenvalue  theorem 
to  impose  the  holonomic  constraints  and  synthesize  the  substructure  and  finally  an 
optimization  technique  to  determine  the  necessary  values  for  substructure  damping  that 
would  fulfill  the  required  bounds  for  system  damping. 
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II.  BACKGROUND 


A.  LAGRANGE 

1.  Theory  and  Derivation 

Lagrange’s  equations  are  commonly  used  to  generate  the  equations  of  motion  for 
second-order  systems  It  utilizes  the  energy  functional  of  kinetic  (T),  and  potential  (U) 
energies,  and  the  damping  dissipation  function  (D)  to  derive  equations  of  motion  for  a 
given  system. 

A  general  fonn  of  Lagrange’s  equation  is, 

d  f  dr)  8U  dD  „ 

- + —  =  QI 

dtydq )  dq  dq 

where  q-q(t)  is  a  generalized  coordinate  of  the  system  and  Q.  is  the  generalized  non¬ 
potential  force  or  moment  resulting  from  excitations  that  adds  energy  to  the  system. 

The  energies  will  be  derived  in  modal  coordinates.  The  following  will  describe 
how  each  components  of  the  Lagrange’s  equation  is  derived  symbolically: 

First  we  must  start  with  the  basic  equation  of  motion  in  physical  coordinates, 

Mx  +  Cx  +  Kx  =  Qj 

where,  M  is  the  mass  matrix,  C  is  the  damping  matrix,  K  is  the  stiffness  matrix  and 
x=x(t).  Using  M  and  K  matrices,  natural  frequencies  ( OJ)  and  mode  shapes  ((j))  were 

calculated  by  solving  the  eigensystem,  |^A-co2M^{(j)}={0}.  In  MATLAB  this 

eigensystem  is  calculated  by  [  ®,  cd\  =  eig(K,  M) ,  where  O  is  a  mass  normalized  modal 
matrix  comprised  of  mode  shapes  cj) .  This  then  results  in  O '  M O  =  M,  =  [I] 
andO'KO=  \K\,  diagonalized  K  matrix  or  a  diagonal  matrix  of  the  structures  natural 
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G-. 

> 

G.. 

1 

frequencies, 

co2 

,  and  finally  ® '  C®  =  [t?]  = 

2^co 

V 

'  V 

,  assuming  that  our 


physical  damping  matrix  is  proportionally  damped.  Since  we  are  currently  working  with 
physical  coordinates,  for  this  thesis  it  is  desired  to  transform  the  equation  of  motion  from 
physical  to  modal  coordinates.  To  proceed  with  the  transformation,  we  need  to  utilize  the 
modal  matrix  ®,  and  the  transfonnation  equation  is  {x}  =  ®{g},  taking  the  first  and 

second  derivatives  with  respect  to  time  t  give  the  following  {x}  =  ®  [cf  and  {x}  =  ®  [if]  . 
As  a  result  the  equation  of  motion  is  now  in  the  form  of, 


r.m  > 

r-.  h 

[I]^  + 

2C,®,7 

q  + 

co; 

V  ‘-V 

v  ‘-J 

Derivation  of  Kinetic  Energy  using  modal  coordinates 

1  1 
Generally  T  =  —  Hq  ,  but  since  M  is  equal  to  [I],  it  can  be  simplified  as  T  =—<7" .  As  a 

result, 


1  N  co 

Z  71  =  1  j= 1 


Therefore,  — 
dt 


rdT ^ 

ySqj 


N  00 

II? 

71=1  7=1 


n(") 

j 


where,  superscript  n  is  the  substructure  designator  and  j  is  the  mode  index 
Derivation  of  Potential  Energy  using  modal  coordinates 
1 

Generally  U =—Kq  Asa  result, 
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^  n= 1  j= 1 

since  [j K]  =  co1 

1  N  co  .  - 

^  n= 1  j= 1 


at/ 


j 

TV  oo 


Therefore,  q{ ^ 

Sq  «=i  j=i 


Derivation  of  Damping  Dissipation  Function  using  modal  coordinates 

1  9 

Generally  D  =  —Cq  Asa  result, 

i  TV  oo  , 

/?  =  1  j= 1 

since  [(?]  =  [2<^Vy] 

i  TV  oo  _ 

"4Z2X!"y"<A 

^  w=l  7=1 

r)D  N  °o 

Therefore,  — 7  =  Z2X!',y",Z"1 

o<Z  n=l  j=l 

Substituting  the  above-derived  energies  into  Lagrange’s  equation  gives  the  following 
result. 


q"]  +  2f"]o.j'l‘iq  "]  -  of”’  q  'n  =0,  since  Qt  =0,  as  a  result  the  jmmode  of  Lagrange’s 


equation  for  the  nth  substructure  is 


--(«)  r,  r{n)  (n)  •  (« )  («)-  (n)  ~ 

q\  ;  +  242  '<»).  ’q\  ’  -  co)  ’  q)  ’  =  0 . 


B.  LAGRANGE  MULTIPLIERS  AND  HOLONOMIC  CONSTRAINTS 


1.  Theory 

Lagrange  Multipliers  (A)  are  included  in  the  generalized  form  of  Lagrange’s 

equation  to  scale  each  holonomic  constraints.  Lagrange  Multipliers’  value  is  a  measure 

of  the  force  of  constraint.  In  terms  of  substructure  coupling,  the  Lagrange  multiplier  is 

interpreted  as  a  measure  of  the  coupling  force  between  substructures.  On  the  other  hand, 
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holonomic  constraints  (/)  ,  are  constraints  on  a  given  coordinate  system;  it  serves  as  an 

interconnecting  condition  between  the  substructures,  and  has  an  overall  effect  of 
diminishing  the  total  number  of  independent  degrees  of  freedom  (N)  to  N  -  r,  where  r  is 
the  number  of  constraints  associated  with  the  interconnecting  region.  Adding  the 
Lagrange  Multipliers  to  the  generalized  Lagrange’s  equation  gives  the  result  as  shown  in 
Klein  [8], 


d_ 

dt 


3U  dD 


(  dT\ 

- - 1 - 1_  _ 

v  dq  )  dq  dq 


e±Kfr 


N  oo 


dq 


=  0  ,  where  fr  =  =  0  and  r=\,...,R  and  the 

n- 1  y=l 


superscript  n  is  the  substructure  designator  and  j  is  the  element  of  mode  shape  or  modal 
coordinate  vector.  This  then  gives  the  final  expression  for  the  equation  of  motion  in  the 
form 


,  R 

••(«)  i  '■>  /'(n)  ( n )  •  («)  («)  ( n )  \  1  /]  n(n)  a 

q)  +  2 Cj  ’0J>)  ]q)  ]  +  <o)  q)  -  2^KPrj  =  0 

r= 1 

Equation  1 

To  better  understand  the  concept  of  interconnecting  region  and  how  it  is  fonnulated  from 
substructures,  a  simple  two  degree  of  freedom  mass  spring  system  as  shown  below  will 
be  used  as  an  example. 


Figure  1.  2  DOF  System  SubsA  (Left),  SubsB  (Right) 

Using  physical  coordinates  to  show  the  connectivity  between  substructure  A  and  B  and 
knowing  that  there  is  only  one  degree  of  freedom  (r  =  1)  that  is  shared  by  the  two 
substructures  after  coupling  the  two;  we  may  now  come  up  with  the  equation  of 
constraint  for  the  two  coordinates  to  be  coupled. 
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fl=x2-x3,  where  x2  and  x,  are  the  degrees  of  freedom  of  each  substructure  that  are 

being  joined  together.  Since  this  thesis  mainly  works  with  modal  coordinates,  we  need  to 
transfonn  each  of  the  degrees  of  freedom  in  terms  of  modal  coordinate,  for  example 

x2  =P<||/<y],l)  +  pg<?g ,  where  P  is  a  component  of  the  modal  matrix  for  each  substructure; 
therefore  the  new  equation  of  interconnecting  region  becomes 


/i  =P!!iVi(1)  +P ilg?  -PuVf2)  -pg#  - 


=  0 


2.  Example  of  Holonomic  Constraint 

Two  objects  each  have  six  coordinates:  x,  y,  z,  r,  0,  (p,  where  x,  y,  and  z  are  the 
center  of  mass  coordinates,  r  is  the  separation  between  the  two  objects,  0  and  (p  are  their 
orientation  in  space.  By  connecting  the  two  objects  with  a  rigid  rod  in  coordinate  r,  we 
are  able  to  decreases  the  degree  of  freedom  both  objects  are  free  to  move  about  by  one; 
therefore,  we  can  say  that  we  placed  the  two  objects  under  a  holonomic  constraint. 

C.  ORTHOGONAL  COMPLEMENT  AND  ZERO-EIGENVALUE 
THEOREM 


1.  Theory 

In  order  to  understand  the  concept  of  orthogonal  complement  we  will  first  define 
the  four  fundamental  subspaces  that  will  be  covered  in  the  explanation  of  orthogonal 
complement.  The  definitions  are  extracted  from  the  book  “Linear  Algebra  and  its 
Applications”  by  Gilbert  Strang. 

Four  Fundamental  Subspaces 

Starting  with  an  m  by  n  matrix  A  we  must  first  define  the  first  part  of  Fundamental 
Theorem  of  Linear  Algebra  Part  1 : 

1 .  If  (At  )  =  row  space  of  A;  dimension  r 

2.  A”  (AJ  =  nullspace  of  A;  dimension  n-r 

3.  \{(Aj  =  column  space  of  A;  dimension  r 
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4.  A”  (AT  )  =  left  nullspace  of  A;  dimension  m-r 

Applying  the  first  part  of  the  fundamental  theorem  of  Linear  Algebra  we  can  now 
talk  about  the  orthogonal  subspace  which  is  a  lead  in  to  the  concept  of  the  orthogonal 
complement.  The  following  are  definitions  from  Strang  [14]  which  talks  about  the 
relationships  between  the  four  subspaces  and  the  relationship  to  the  theory  of  the 
orthogonal  complement. 

1.  Two  subspaces  V  and  W  of  the  same  space  R"  are  called  orthogonal  if 
every  vector  v  in  V  is  orthogonal  to  every  w  in  W  :  vTw  =  0  for  all  v  and 
w. 

2.  For  any  m  by  n  matrix  A,  the  nullspace  A”  (A)  and  the  row  space  If  ( A 1 ) 
are  the  orthogonal  subspaces  of  R" .  Similarly,  the  left  A”  ( At )  and  the 
column  space  1 1(A)  are  orthogonal  subspaces  of  R"' . 

3.  Orthogonal  complement  definition:  Given  a  subspace  V  of  R" ,  the  space 
of  all  vectors  orthogonal  to  V  is  called  the  orthogonal  complement  of  V, 
and  denoted  by  V1 . 

Using  this  terminology,  the  nullspace  A '(A)  is  the  orthogonal  complement  of 
If  ( AT ):  < A  (A)  =  ( If  (At)  )±.  At  the  same  time,  the  opposite  relation  also  holds:  The 
row  space  If  ( AT  )  contains  all  vectors  that  are  orthogonal  to  the  nullspace.  This  is  not  so 
obvious  from  the  construction,  since  in  solving  Ax  =  0  we  started  with  the  row  space  and 
found  all  x  that  were  orthogonal  to  it;  now  we  are  going  in  the  opposite  direction. 
Suppose,  however,  that  some  vector  z  in  R"  is  orthogonal  to  the  nullspace  but  is  outside 
the  row  space.  Then  adding  z  as  an  extra  row  and  A  would  enlarge  the  row  space 
without  changing  the  null  space.  But  we  know  that  there  is  a  fixed  formula: 
r+(n-r)=n  ,  or  dim(row  space)  +  dim(nullspace)  =  number  of  columns.  Since  the  last 
two  numbers  are  unchanged  when  the  new  row  z  is  added,  it  is  impossible  for  the  first 
one  to  change  either.  We  conclude  that  every  vector  orthogonal  to  the  nullspace  is 

already  in  the  row  space:  ^  ^  )  =  A” (A)1  (Strang)  [14]. 
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The  same  reasoning  applied  AT  produces  the  dual  result:  The  left  nullspace 
eN”  ( A 1 )  and  the  column  space  R  (A)  are  orthogonal  complements  of  one  another  in  W" . 
This  completes  the  second  half  of  the  fundamental  theorem  of  linear  algebra.  The  first 
half  gave  the  dimensions  of  the  four  subspaces,  including  the  fact  that  row  rank  =  column 
rank,  and  now  we  know  that  they  are  not  only  perpendicular;  they  are  orthogonal 
complements  (Strang)  [14]. 

Zero-eigenvalue  theorem  (Walton  and  Steeves)  [16]:  Let  E  be  n  x  n  matrix  AT  A  _ 
If  the  rank  of  E  is  m,  therefore  the  eigenvalues  of  E  is  also  m  and  it  has  s  =  n  -  m  “zero 
eigenvalues”.  These  zero-eigenvalues  have  corresponding  eigenvectors,  which  is  called 
the  nullspace  of  matrix  E.  Having  solved  for  the  set  of  eigenvectors  corresponding  to  the 
zero-eigenvalues  of  matrix  E  we  can  call  these  set  of  eigenvectors  T,  which  is  an  n  x  s 

matrix.  Moreover,  by  post  multiplying  the  equality  E  =  Ar  A  wjth  the  newly  generated  T 
matrix  we  will  come  up  with  the  expression  E  T  =  AT  A  T  =  0  and  by  premultiplying  the 
said  equation  with  Tt  ,  the  equality  remained  equal  to  zero  as  shown  by  Tr  AT A  T  =  0; 
further  matrix  manipulation  will  simplify  this  term  to  (AT)7  AT  =  0,  which  finally  leads 
to  the  fact  that  matrix  T  is  an  orthogonal  complement  of  matrix  A  (Kamman)  [7], 

The  two  theorems  discussed  above  is  one  of  the  central  aspects  that  will  be 
employed  in  this  thesis;  the  above  said  matrix  T  ,  will  be  utilized  as  a  transfonnation 
matrix  of  which  will  enable  us  to  generate  coupled  system  equations  of  motion  given 
only  modal  properties  of  each  substructure. 

2.  Example 

MATLAB  code  was  generated  to  show  the  usefulness  of  zero-eigenvalue  theorem 
in  substructure  coupling  using  constraints.  Using  the  example  from  William  Walton  and 
Earl  Steeves’  technical  report,  the  MATLAB  code  is  as  follows: 

clear 

clc 

%Mass  Spring  System 

%M  1  ==/V==M2==/V==M3==/V==M4==/V==M5 

%Example  from  Walton  and  Steeves  paper 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
ml  =  [1,1, 1,1,1];  %mass  values  for  subs  1 

kl  =  [1,1, 1,1];  %stiffness  values  for  subs  1 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

K=zeros(5,5); 

%Building  of  K  and  M  matrices  for  each  substructure; 

index  =  [0 , 1]; 
for  i  =  1  :  length(kl); 
index  =  index  +  1 ; 

ka  =  [k  1(1  ,i)  ,  -k  1(1  ,i)  ;  -kl(l,i)  ,  kl(l,i)]; 

K(  index,  index)  =  K(index,index)  +  ka; 
end 

M  =  diag(ml); 

[Phi,Lam]  =  eig(K,M); 

W  =  sqrt(Lam) 

%Subdivided 
Kbar=zeros(8,8); 
index  =  [-1  ,0]; 
for  i  =  1  :  length(kl); 
index  =  index  +  2; 

ka  =  [kl(l,i) ,  -kl(l,i) ;  -kl(l,i) ,  kl(l,i)]; 

Kbar(index,index)  =  Kbar(  index,  index)  +  ka; 
end 
Kbar; 

mlbar  =  [1,. 5, .5, .5, .5, .5, .5,1]; 

Mbar  =  diag(mlbar); 

%Equation  of  Constraint 
C  =  [0, 1,-1, 0,0, 0,0,0; ... 

0,0, 0,1, -1,0, 0,0; ... 

0,0, 0,0, 0,1, -1,0]; 

E  =  C'*C; 

[U,lam]  =  eig(E); 

T  =  null(E); 

%System  M  and  Stiffness  matrices 
MBAR  =  T'*Mbar*T 
KBAR  =  T'*Kbar*T 

[Phibar,Lambar]  =  eig(KBAR,MBAR); 

Wbar  =  sqrt(Lambar) 
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Wbsr  = 


0  +  0.  OOOOi  0 

0  0.  6180 
0  0 

0  0 
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0 

0 

1.  1756 
0 
0 


0 

0 

0 

1.  6180 
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0 

0 

0 

0 

1.  9021 


The  natural  frequencies  W  (solution  using  equation  of  motion  for  the  system)  and  Wbar 
(solution  from  subdividing  the  system  into  5  parts)  matches. 

D.  DECOUPLING  APPROXIMATION 


Most  structures  are  analyzed  via  classical  damping  models,  where  physical 
damping  matrix  C  is  diagonalized  by  calculating  the  triple  product  O'  CO  =  [2(^0);] 

which  enables  us  to  extract  modal  damping  ratios  for  the  structure  being  analyzed.  On 
the  other  hand,  some  structures  have  a  dynamic  characteristic  that  is  non-classically 
damped,  meaning  they  are  not  diagonalized  by  the  mode  shapes  of  the  structure;  as  a 
result,  extracting  the  modal  properties  of  the  structure  becomes  a  tedious  task.  Structures 
that  exhibit  this  kind  of  characteristics  are  those  that  come  from  the  coupling  of 
substructures,  structures  that  have  different  damping  characteristics  or  structures  that  has 
highly  damped  components  (Xu)  [17].  Fortunately,  we  may  circumvent  the  problem  of 
having  a  non-classically  damped  structure  by  applying  the  concept  of  decoupling 
approximation,  which  allows  for  the  decoupling  of  structures  modal  matrix  by  simply 
ignoring  the  off-diagonal  components  of  the  matrix  and  using  the  main  diagonal  elements 
to  come  up  with  the  damping  ratios  of  the  substructure.  In  order  for  us  to  implement  the 
decoupling  approximation  for  the  structure  in  hand,  we  must  first  satisfy  some  necessary 
requirements,  the  following  will  explain  how  we  formulate  required  values  such  as 
second-order  eigenvalues  and  eigenvectors  and  see  if  they  comply  with  the  necessary 
requirements. 

Starting  with  the  general  equation  of  motion  Mox  +  Cax  +  Kox  =  0,  as  previously 

shown  above  we  are  able  to  generate  a  mass  normalized  modal  matrix  ®o  by  solving  the 
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eigensystem  ^K0 -a>2joMo ]{®o}={0}  in  MATLAB.  Given  ®o,  and  the  relationship 

x=®0<7,  we  can  now  convert  the  equation  of  motion  in  modal  coordinates.  Therefore 
the  resulting  equation  is  Iq  +  Gq  +  A q  =  0 ,  where  /  is  an  identity  matrix,  G  is  the  modal 
damping  matrix,  which  was  not  diagonalized,  and  A  is  the  decoupled  stiffness  matrix. 
The  next  step  is  to  solve  for  the  second-order  eigenvalue  problem  corresponding  to  the 
modal  equation  of  motion  derived  above.  The  formulation  of  the  second-order 
eigenvalue  problem  is  shown  in  the  equation  [x:  I  +  Sj  C  +  A]HA  =0 ,  where  ,s';  and  HA 

are  the  system’s  eigenvalues  and  eigenvectors.  These  HA ,  of  which  are  mass 
normalized  modal  matrix,  are  then  used  to  obtain  the  structure’s  mode  shapes  ®  in  the 

original  displacement  coordinates  x  by  the  relation  ®y  =  ^  i//vO m  ,  where  y/  is  the  zth 


component  of  HA .  Moreover,  the  generated  eigenvalues  s  .  are  complex  natural 
frequencies  that  are  comprised  of  structures  natural  frequencies  coj  and  f,  and  can  be 

expressed  in  thee  form  sf  =  oy  |  -  C;  ±  i^l  -  j .  To  find  out  whether  or  not  the 

structure  that  is  currently  being  worked  can  be  accurately  represented  by  the  decoupling 
approximation,  we  must  see  if  the  generated  mode  shapes  HA  and  eigenvalues  s  . 

satisfies  the  requirement  y//  =0(s)  for  all  values  of  z  ^  j  and 


s  =  (0 

JO  JO 


^z±iV1-^)’where  = 


e„ 


{<°io+G>Jo) 


and  G.  is  the  zyth  component  of 


the  damping  matrix  G  and  0(6')  represents  a  tenn  of  the  same  order  of  magnitude  as 
s  or  smaller,  and  s  is  relatively  smaller  than  one  (Xu)  [17]. 
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E. 


GERSCHGORIN’S  THEOREMS 


1.  Theory 


The  Gerschgorin  disk  theorem  is  a  theorem  which  is  used  to  bound  eigenvalues  of 
a  square  matrix.  This  was  shown  in  the  book  by  Leonard  Meirovitch  [12]  in  a  more 
complete  detail,  the  proof  is  as  follows: 


“Let  us  write  the  eigenvalue  problem  associated  with  the  n  x  n  matrix  A  in  the 

n 

form  'y\akjxj  =  Xxk,  k  =  l,2,...,n  next,  let  us  assume  that  xm  is  a  component  of  the  vector 

j= i 

x  with  the  largest  modulus,  |xm|=max|jcy|  ( /  =  1,2,...,/?) ,  so  that  letting  k  =  m  in  above 


equation,  we  can  write  (X  -cinm)xm  =  'Yjamjxj  Hence, 

j= i 

j^m 


Dividing  the  inequality  through  by  \xm\,  we 

i= i  j=i 

jiXm  j^m 


obtain  L  -  amm  <  ^  \amj  .  This  inequality  is  called  the  first  theorem  of  Gerschgorin  and 

j= i 

j^m 


can  be  stated  as  follows:  Every  eigenvalue  of  the  matrix  A  lies  in  at  least  one  of  the 


circular  disks  with  centers  at  a.„m  and  radii  r  =  V  a  .  .  The  disks  are  sometimes 

mm  m  mj 

j= i 

j^m 


referred  to  as  Gerschgorin’s  disks.”  Meirovitch  [12].  This  theorem  seemed  useful  in 
coming  up  with  possible  bounds  for  both  system  and  substructure  damping  but, 
unfortunately  for  this  theorem  to  be  of  much  use  in  estimating  modal  damping,  it  is 
important  that  the  off-diagonal  elements  of  matrix  A  be  small  relative  to  the  main 
diagonals,  as  a  result  this  theorem  was  not  put  to  use  in  this  thesis. 
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2.  Example 


From  Meirovitch’s  book 

"2.5  -1  0  ' 

Consider  the  matrix^  =  -1  5  -  sfl  ,  the  centers  of  Gerschgorin’s  disk  are  the 

0  -V2  10 

diagonals  an  =2.5,  a22  =5,  a33=10  and  the  corresponding  radii  are 

i\  =  |a12|  +  |a13|  =|-l|  + 10|  =1 

f2  =  |^?21  |  1^23  |  =|_  l|  =  1  +  V2  , 

r3  =|a31 1  +  |fl32  |  =|0|  + 

The  eigenvalues  for  matrix^,  areXj  =2.1 19322,  X2  =5,  X3  =10.380678  . 


Figure  2.  Gerschgorin’s  Disk  shows  how  each  eigenvalues  falls  within  each  disk 

F.  METHODS  FOR  PRODUCING  PROPORTIONALLY  DAMPED 
SUBSTRUCTURES 

1.  Method  A 

Modal  damping  ratio  (Q  is  a  physical  property  of  a  system  that  is  typically 
measured  in  a  vibration  test  experiment.  When  conducting  structural  coupling  we  can 
assume  value(s)  for  a  given  structure’s  C,  for  each  corresponding  modes,  a  typical 
minimum  value  used  commonly  used  is  around  two  percent  (C  =  0.02).  In  doing  so,  we 
assume  the  given  structure  to  be  a  proportionally  damped  system. 
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2. 


Method  B 


Typically  the  way  we  would  come  up  with  proportional  damping  is  by  making  the 
physical  damping  matrix  proportional  to  both  mass  (M)  and  stiffness  (K)  matrices  with  a 
given  proportionality  constants  a  and  [3.  C  =  aK  +  [}  M  ,  with  this  proportionally 
damped  C  matrix  we  may  now  extract  damping  ratios  by  first  diagonalizing  the  physical 
damping  matrix  C  by  conducting  the  triple  product  O'  CO  =  [2<^icoi  ] ,  where  O  is  the 


modal  damping  generated  from  solving  the  eigensystem  consisted  of  K  and  M  matrices 
and  by  dividing  the  corresponding  values  of  the  diagonal  matrix  with  the  corresponding 


natural  frequency  multiplied  by  two  in  the  form,  C,  = 


2C,co, 

2co(. 


Although  we  may  create  proportionally  damped  substructures  by  providing  values 
for  C,  or  extracting  C  from  a  proportionally  damped  physical  C  matrix,  it  is  not 
guaranteed  that  the  final  damping  matrix  C  will  be  a  proportionally  damped  one;  in  fact 
coupled  system  in  general  are  those  of  the  non-proportional  type  Udwadia,  F.  E.  [15]. 
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III.  MILESTONES  OF  THESIS  WORK 


A.  COUPLING  OF  MASS  SPRING  SYSTEM 

Synthesizing  a  simple  mass  spring  system  will  be  used  to  illustrate  the  process  of 
constructing  system  matrix  and  extracting  from  it  system  natural  frequencies  and  modal 
damping.  Figure  1  below  is  the  model  that  is  being  synthesized. 


Figure  3.  2  DOF  System  SubsA  (Left),  SubsB  (Right) 

The  following  will  explain  the  three  different  methods  that  were  used  to  calculate 
system  physical  properties.  Furthermore,  coding  to  come  up  with  physical  property 
values  was  done  in  MATLAB  and  since  the  method  of  constructing  substructure  A 
(subsA)  and  substructure  B  (subsB)  matrices  are  identical,  this  thesis  will  focus  only  on 
substructure  A. 


1.  Method  1;  Basic  Structural  Coupling 


cl  c2  c3  c4 

Figure  4.  3  DOF  Mass  Spring  System 
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Subs  A 


Initially,  physical  properties  Mass  (M),  Stiffness  (K)  and  modal  damping  ratio 
(£)  for  each  mode  are  given.  Using  M  and  K  matrices,  natural  frequencies  {(»)  and 

mode  shape  (®)  were  calculated  by  solving  the  eigensystem,  |^if-co2Mj{(j)}={0}.  In 
MATLAB  this  eigensystem  is  calculated  by  [  ®,  co\  =  eig(K,  M) ,  where  O  is  a  mass 
nonnalized  modal  matrix.  This  then  results  in  ®'M®  =  [I]  and®’K®  =  [/f ] , 

diagonalized  K  matrix.  Since  co  is  now  known,  physical  damping  matrix  (C)  can  now  be 
solved  via  the  equality  ® 'C®  =  [2£co] ,  where  [2Cco\  is  a  diagonal  matrix,  this  will  be  a 
proportionally  damped  substructure  since  modal  damping  is  initially  given.  Finally 
solving  for  C,  since  ®'M®  =  [I]  and  ®'C®  =  [2^co]  we  solve  for  the  physical  damping 
by  conducting  the  following  matrix  manipulations, 

C  =  ®  [2<^co]®  ’ 
since  M®=®  [I] ;  M®  =  ® 
and  ®'M=  [I]®’;  ®’M=  ®’ 
therefore, 

C  =  M  ®[2<^co]®'M 


The  mass,  stiffness  and  damping  matrices,  M,  C,  and  K,  for  substructure  B  are 
calculated  in  the  same  way.  Having  created  Ma,  Ka,  and  Ca  for  SubsA  and  Mb,  Kb,  and 
Cb  for  SubsB,  the  systems  Ms,  Ks  and  Cs  can  be  calculated  by  adding  the  terms  of  the 
substructure  matrix  that  corresponds  to  the  same  position  in  the  system  matrix.  For 


example, 


Mau 

Ma21 

0 


Man  0 

Ma22+Mbn  Mbn 
Mb2]  Mb22 


Having  calculated  for  system’s  physical  mass,  stiffness  and  damping  matrices  the 
next  step  is  to  solve  the  eigensystem  [Ks -<zrMs  jj®sj=jO} ,  where  the  subscript  (s) 

denotes  that  we  are  making  use  of  systems  mass  and  stiffness  matrices,  to  detennine  the 
values  for  cos  and  ®s.  With  these  two  modal  properties  in  hand  we  are  now  be  able  to 

solve  for  system  modal  damping  ratios;  in  order  to  do  so  the  triple  product  ®S'CS®S  is 
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conducted  in  order  to  modalize  the  damping  matrix  which  then  enables  us  to  extract  said 
damping  ratios  from  the  matrix  Os 'Cs  Os  =  [2^scoJ  .  As  mentioned  above  in  the 
decoupling  approximation,  we  will  ignore  the  off-diagonal  terms  of  the  modal  mass 
matrix  when  calculating  for  the  damping  ratio  C,s ;  this  assumption  was  made  because  we 
are  not  guaranteed  to  have  a  proportionally  damped  system  via  coupling  of  two  sub 
substructures.  (Xu)  [17].  Therefore,  =  C;i  / 2coi; ,  where  i  =  1,2,3,  to  N  and  ii  is  the 
diagonal  value  for  each  corresponding  system  matrix. 

2.  Method  2;  Calculating  System’s  Damping  Ratio  and  Natural 
Frequencies  by  Transformation  Method  Using  Zero-Eigenvalue 
Theorem 


In  this  section  structural  synthesis  is  conducted  by  implementing  the  Zero- 
Eigenvalue  Theorem.  This  theorem  will  provide  the  necessary  transfonnation  matrix  that 
will  convert  a  block  matrix  of  substructures  A  and  B  to  a  synthesized  system  structure. 

Using  Lagrange’s  equation,  Lagrange  multiplier  and  holonomic  constraints  the 
final  equation  of  motion  is  in  the  form  i/'f1  +  2^{"] (o"] q"]  +  co" 1  q"]  -  ^ Z  /f " 1  =0  for 

r= 1 

each  substructure. 

This  example  will  illustrate  how  method  two  is  applied  to  the  above  mentioned 
mass  spring  system.  Using  the  generalized  equation 

'qf  +  -2X4"’  =  0 


Each  substructures  equations  of  motion  are  as  follows 


Substructure(A) 
q\l)  +  2  +  ©}’ 

q^  +  2£^a$c$  +  cOj 
Substructure^ ) 


-VS)  =  o 

-4A(J=0 


#+2  ?r+w= 0 

+  2  £M2)£]  +  <1?  +  M(2}  =  0 


'(2)  ,  ^(2)  A1) 


(2) 
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The  next  step  is  to  put  the  above  equations  of  motions  in  the  form 

q  +  ICoxj  +  co2q-[A]T  A  =0 ,  where  A  =  -/?1,l]i  -  /?' /3\]]  ,  a  matrix  of  holonomic 

constraints.  The  block  diagonal  consisted  of  2 C,a>  shall  be  called  blockC  and  co2  as 
blockK  for  labeling  purposes. 

1  0  0  0  1U«1  \ap  0  0  0  lL«l  ' 

m  0  2CiV1)  0  0  qf  0  0  0  qf  -0$ 

WU  +  o  o  o  #  r  „  „  4*0  Ur  $  w 

M  [  0  0  0  [o  0  0  cof  JU2)J  U? 

Now  let  D  =  ArA ,  and  then  find  the  orthogonal  complement  of  D  matrix  by  using  the 
“null”  function  in  MATLAB,  which  gives  the  basis  vectors  for  the  null  space  of  D 
matrix.  These  basis  vectors  will  now  be  called  matrix  N,  this  N  matrix  will  be  used  as  a 
transfonnation  matrix  to  transform  the  block  diagonal  matrices  blockC  and  blockK  to 
system  modal  damping  matrix  and  modal  stiffness  matrix.  After  the  transforming  the 
block  matrices  to  system  matrices  the  systems  natural  frequencies  and  damping  ratios 
were  then  extracted  from  the  system’s  modal  damping  and  modal  stiffness  matrices.  The 
following  will  show  the  process  on  how  to  conduct  the  transformation. 

First,  as  discussed  above,  we  must  extract  the  block  matrices  of  oJ"'1  and 

oJj1'  from  the  above  equations  of  motion;  the  upper  block  corresponds  to  SubsA  and  the 
lower  block  to  that  of  SubsB 

24'1(V1)  0  0  0 

o  2  0  o  =  [blockC]  4x4 

0  0  2 0 

0  0  0  2 4fW2) 
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0  0  0 


0  cof  0  0 

0  0  cof  0 

0  0  0  cof 


=  [ blockK ]  4x4 


Next,  form  a  triple  product  to  transform  matrices  blockC  and  blockK  using  matrix  N  to 
generate  the  new  matrices  G  and  L  . 

IG]  =  Nt  [blockC]  N  3x3 

=  NT[b!ockK]N  3x3 

With  the  above  given  matrices,  we  can  now  put  it  in  the  first  order  form  by 
transfonning  the  second  order  differential  equation  to  first  order 


[«;■»] 

w\ 

"o  r 

differential  < 

••( n ) 

Wj  j 

>  =  [£]< 

•  («) 
W  J 

> ,  where  E  is  equal  to  the  block  matrix,  E  = 

L  G 

next  step  is  to  solve  for  the  complex  eigenvalues  of  matrix  E,  ^i=oi  ±  Qd  j ,  where 
a;  =  <^03;  and  Qd  is  the  undamped  natural  frequency  of  the  system.  Finally,  using  the 

relationships,  Od  =  co^/l  -  Q*  and  a;  =  C,^,  we  can  now  solve  for  systems  modal 
damping. 

The  resulting  co  and  Q  from  this  method  matches  the  results  coo  and  of 
method  one  above. 


3.  Method  3  Branches  Out  of  Method  2 

Another  method  that  was  conducted  to  get  values  for  co  and  tT  branches  out  from 
that  of  the  transfonnation  method.  In  this  method,  we  make  use  of  the  already  generated 
[G]  and  [L]  matrices.  First,  system  natural  frequencies  co  and  mode  shapes  ®  are 

extracted  from  [F]  by  solving  for  the  eigenvalues  and  corresponding  eigenvectors.  Upon 

extraction  of  [F]’s  eigenvectors,  the  next  step  is  to  take  the  triple  product 
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New_G  =  0'[G]0  .  The  matrix  New_G  will  now  serve  as  the  modal  damping  matrix 

for  the  system,  though  it  is  not  a  diagonal  matrix,  since  it  is  not  proportionally  damped, 
we  are  once  again  making  use  of  the  decoupling  approximation  to  calculate  the  values 
for  the  system  modal  damping.  The  final  step  of  this  method  is  to  individually  extract  the 
diagonal  values  of  New_G  matrix  and  divide  it  with  the  corresponding  system’s  2co , 


c,= 


2<>, 

2  co. 


by  doing  so  we  can  make  an  observation  that  the  solution  will  indeed  be 


equal  to  that  of  the  above  calculated  values  from  the  first  two  methods. 


4.  Comparing  the  Calculation  Cost  Between  Method  2  and  3 

Since  computational  cost  is  a  big  issue  in  solving  for  numerical  problems.  We 
need  to  investigate  how  to  improve  computation  time  for  almost  any  task  that  is  being 
done  computationally.  As  a  result,  for  this  thesis  we  need  to  look  at  and  compare  which 
method,  method  two  or  three,  would  generate  a  solution  in  a  more  economical  way.  In 
the  two  methods  discussed  above,  the  main  difference  in  the  methodology  was  the  size  of 
the  matrix  that  is  being  solved.  Method  two  used  a  matrix  of  dimension  2 n  by  2 n  in 
solving  for  the  solution  while  method  three  made  use  of  a  matrix  of  size  n  by  n.  Knowing 
this  fact,  we  further  investigated  which  matrix  takes  longer  to  calculate.  The  first  part  of 
the  investigation  dealt  with  calculating  for  the  time  it  takes  to  compute  for  the 
eigenvalues  symmetric  matrices  of  varying  sizes  starting  from  size  n  by  n  to  16 n  by  1 6n. 
Using  a  generated  MATLAB  code  created  by  Gordis,  it  is  clear  that  the  larger  the 
dimension  of  a  matrix  becomes  the  more  time  it  takes  to  calculate  for  the  eigensolution. 
The  second  part  of  the  investigation  called  for  finding  the  exact  relationship  between  the 
cost  (time)  and  the  sizes  of  matrices  that  are  being  solved  for  its  eigensolution.  Since  we 
are  employing  MATLAB  to  solve  for  the  eigenvalue  problem,  which  makes  use  of  QR 
and  QL  algorithms,  where  Q  is  the  orthogonal,  R  is  the  upper  triangle  and  L  is  the  lower 
triangle  of  any  matrix  A;  we  can  compute  the  workload  of  a  general  matrix  to  be 
approximately  0(n  ).  This  then  makes  method  three  a  better  choice  in  solving  for  the 
solution. 
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B.  INITIAL  WORK  AND  THE  REASON  FOR  IMPLEMENTING  ZERO 
EIGENVALUE  THEOREM 

The  main  motivation  on  why  this  concept  is  being  implemented  on  this  thesis  was 
the  fact  that  it  made  synthesis  of  two  substructures  less  tedious  in  an  algebraic  stand 
point.  During  the  initial  phase  of  the  thesis  work,  the  method  of  choice  in  synthesizing 
the  two  substructures  was  to  transform  the  equations  of  motion  of  each  substructure  from 

a  second  order  differential  equation  into  a  first  order  differential,  \  f  =  [dll  ^  I  where 

wJ  m 

[A]  is  the  system  matrix  which  will  be  used  to  calculate  for  system’s  dynamic  response  . 

To  generate  this  system  matrix  we  must  first  consider  the  equations  of  motion  from  each 
substructure  that  was  generated  via  Lagrange’s  equations.  Using  the  same  mass  spring 
system  that  was  mentioned  above  and  its  equation  of  motions 

q[l)  +  2£(%(1)g1W  +  q\l)  - =  0 

q(i  +  2£lVV1)  +  4°  q(i  ~  =  0 

q[2)  +2 Cl2)co[2)q[2)  +oj{2]  q{2)  +  A^2)  =0 
q(2)  +2 Ci2)co{2]q{2)  +  a>(2)  q(2)  +  =  0 

then  solving  for  <y(;")  gives  the  following  equations 

9,0)  =  ~  (EOM\) 

qf=-2  +  AA(S  (EOM2) 

if  -V!2J  (EOM4) 

Using  the  interconnecting  conditions  between  the  two  substructure  components 

N  oo 

V/'*  =  0  and  r  =  l,...,R  ,  since  r  =  1,  the  explicit  form  is 


The  next  step  is  to  take  the  first  and  second  derivatives  of  the  constraint  equation  with 
respect  to  the  modal  coordinate  q  which  results  in  the  following  equations 


A  =  fiSq?  +  /$$  ~  Aff#  -  PxW?  =  0 
A  =  AW}  +  PS^  ~  P$q?  ~  A<22^22)  =  0 
A  =  fiS'q?  +  PS&  -  Pu'4  !2)  -  Pi- Td?  = 0 


Then  substitute  in  the  corresponding  values  of  cji")  from  the  equations  of  motion  above 
and  then  solve  for  the  Lagrange  multiplier  A ,  stepping  through  the  process: 

PSqP+pQq?  -Pu’q?  -PiTit?  =0 

=>  Substitute 


A(,?  (~2rfI)^(1)4ri(1)  -^i(I)  ?i(1)  +^)+P$[-^)4)$-4)  <l2)  +  \Pi} 

-fit?  t-tfVW  -  <TV!)  -  k?  Md  V: -  d2)‘??> 1  -  2,  A 


(2)1  _ 


=  0 


=>  Distribute 

+  A?2dV2)«i2) +-W42’ MAf/  =0 


Solve{/ 1) 


+a?  +a,(2) 


1,1 


1,2 


♦aS^P  -  Aa)2rfa  w  -a<M2)2# 


a 


aSW'V 


A^+A^A^f 


At?  +A?/  +Ai?  +Aj? 

ASrtMW 


u  +  Au  +  Au  +  Au 


Af>f)2# 


a/?2  +aS2  +a,(;)2  +Af/ 


A5WV 


Au)2  +  Ai(1/  +  Au)2  +  A/2/ 


+A?/)" 


Au  +  Au  +  Aii  +  Au 
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"(2) 

9i  = ' 


aS'aS’2  dV1 

A/(1)  - 

aMM1’’ 

(Af +A"/+A(f +A1?) 

Cl\ 

(Ad+A'C'+Af+A12/) 

(i) 


AfAlS  2?M] 


S  +k'.  +A.r  +kJ 
|,)’  +  a!1/  + aS|! + a(2): 


ii’ 


1,1 


1,2 


«r2)- 


1,1 


1,2 


-2cl2)42) 


a\2)  + 


Pit 


-  CO, 


,(2)2 


9i 


(2) 


(if 

1,1 


+  A2’+A?’+A,2> 


1,2 


A'j’A'M1' 


»  +A('/+aJ?,+A<2)' 


U 


'1,2 


^2 


(2) 


"(2) 

^2  = 


AdACdd0®,"1 

aM_ 

AW 

> 

tf+Af+Af+A?) 

Cl\ 

(Aj1’  +A?  +a(?’  +a<2/) 

9i 


(1) 


tf+AS+Aif+Af 


ii’ 


aIMM" 


<*+/&+/&+/& 


1,1 


1,2 


’u  +A',2  +A',iJ  +A 


?w 


(2f 


?(2)2 

1,2 


-2rf  ®; 


<2K.f> 


</!2)  + 


(?  +A',2  +AV  +  A 


?w 


(2f 


f2)2 

1,2 


-  61 


,(2f 


?i 


(2) 


/gdgfag 

"l+Af+A'f  +  A(2)’ 


1,1 


1,2 


-2d2’®;2’ 


«f  + 


a|2/«42)‘ 


ff  +C  +A.;’  +A',: 


)W 


(2f 


?<2f 


1,2 


-<Z>- 


,(2)2 


(2) 

?2 
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Finally,  assemble  the  system  matrix  .4, 


[?!"] 

K’l 

•  (2) 

«,(3) 

•  (2) 

•  =  M- 

(2) 

q\ 

?!" 

Ai] . 

As  we  can  see  with  only  one  interconnecting  region  the  calculation  process  takes  an 
exorbitant  amount  of  time.  As  a  result  the  process  of  using  zero-eigenvalue  proved  to  be 
a  highly  beneficial  method  of  synthesizing  substructures. 

C.  CONSTRAINT  ON  THE  SYSTEM  MODAL  DAMPING 

The  method  of  substructure  coupling  was  put  in  place  in  order  to  determine  what 
some  of  the  possible  bounds  are  for  system  damping.  With  this  concept  in  mind,  we  may 
now  focus  on  the  aspect  of  predicting  the  possible  bounds  for  sub  structural  damping  that 
would  suffice  the  necessary  damping  ratio  for  the  system. 

Knowing  what  the  desired  range  of  system  damping  should  be,  of  which  we  can 
refer  to  as  Zetamax  and  Zetamin,  and  having  calculated  systems  response  from 
substructure  mass  and  stiffness  matrices  from  either  methods  two  or  three  mentioned 
above;  we  may  now  proceed  in  calculating  for  the  main  motivation  of  this  thesis  which 
are  the  necessary  bounds  for  sub  structural  damping  that  would  not  violate  the  necessary 
bounds  for  system  damping.  This  task  can  be  achieved  by  slightly  modifying  the  above 
methodology  of  substructure  coupling.  The  whole  process  is  explained  below: 

1.  First,  by  choosing  either  method  two  or  three  above,  we  must 
calculate  substructure  response,  given  initial  values  for  substructure 
mass  and  stiffness  matrices.  As  mentioned  before,  by  solving  the 
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eigenvalue  problem  [if-co2M]{(j)}={0},  we  will  come  up  with  the 
substructure’s  natural  frequencies  ( OJ)  and  mode  shapes!  ®  ). 

2.  Given  substructure  mode  shapes,  Oa  and  ®B  we  are  now  able  to 
extract  from  these  mode  shapes  the  values  for(p) ,  which  are  elements 

of  each  mode  shape  that  corresponds  to  each  interconnecting  points 
for  coupling  the  substructures,  to  help  fonnulate  the  constraint 
matrix  [A] . 

3.  This  constraint  matrix  will  then  be  utilized  in  forming  the 
transformation  matrix  using  zero-eigenvalue  theorem.  This  is  shown 
by  first  solving  for  D  =  ATA,  then  by  taking  the  orthogonal 
complement  of  D,  we  can  generate  matrix  N  which  is  consisted  of  the 
basis  vectors  for  the  null  space  of  matrix  D.  This  will  then  allow  for 
the  calculation  of  system’s  natural  frequencies  via  transfonnation  of 
blockK,  which  is  a  block  diagonal  matrix  consisted  of  substructure’s 
natural  frequencies  squared  ( or  j .  The  output  from  this 

transformation  is  exactly  what  matrix  L  is  from  method  three  above, 
with  this  matrix  we  can  come  up  with  systems  natural  frequencies  and 
mode  shapes  by  solving  for  the  eigenvalues  and  eigenvectors  of 
matrix  L  in  MATLAB. 

4.  After  generating  system  natural  frequencies  cos  and  mode  shapes  ®s, 
and  prescribing  values  for  system  modal  damping  which  falls 
between  desired  values  Zetamax  and  Zetamin,  we  can  come  up 
with  system’s  physical  damping  matrix,  via  the  equation 
C,  =  ®,[2£>,]®/ 

5.  Given  physical  damping  matrix  Cs  and  the  previously  calculated 
transformation  matrix  N  we  can  decompose  this  system  matrix  into 
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two  matrices  consisted  of  block  diagonal  blockC  consisted  of  2 C,a> 
for  each  substructure. 


6.  Finally,  since  the  natural  frequencies  and  mode  shapes  were 
calculated  early  on  for  each  substructure,  extraction  of  each 


substructures  modal  damping  can  be  perfonned  via  = 


2<>, 

2co. 


.  This 


then  enables  us  to  determine  the  bounds  for  each  substructure 
damping  that  would  satisfy  the  desired  Zeta  min  and  Zetamax. 


D. 


OPTIMIZATION 


1.  Work 

Now  that  we  have  an  idea  what  the  necessary  bounds  are  for  substructure 
damping  given  a  constraint  on  system  damping.  In  this  portion  of  the  thesis  we  are 
tasked  to  implement  the  method  of  optimization  to  ease  the  process  in  calculating  for  the 
values  for  sub  structural  damping,  bound  by  a  lower  and  upper  limit,  which  would  satisfy 
the  given  constraints  for  system  damping  for  all  modes. 

Instead  of  formulating  our  own  code,  a  built  in  MATLAB  function  FMINCON  is 
put  to  use  in  finding  these  substructure  damping.  In  MATLAB,  fmincon  is  a  built  in 
function  which  aids  in  finding  a  constrained  minimum  of  a  given  objective  function  of 
several  design  variables  starting  at  an  initial  estimate.  This  is  generally  referred  to  as 
constrained  nonlinear  optimization  or  nonlinear  programming. 

Fmincon  will  minimize  the  objective  function  which  is  a  summation  of  all 

N 

substructure  damping  rations,  the  design  variables  (DV),  f(C)  =  ZX  ,  subject  to  a  lower 

i= 1 

and  upper  bound  for  the  solution  (DV)  and  a  nonlinear  constraint  equation  which  then 
provides  upper  and  lower  bounds  for  system  damping. 
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E. 


EXAMPLES  OF  WORK  THAT  WAS  DONE 


1.  Mass  Spring  System 

All  of  the  work  done  with  this  type  of  system  was  discussed  above. 

2.  Beam  Structure 


Subs  A 


SubsB 

Figure  5.  Cantilever  Beam 


7 


Figure  6.  Synthesized  Cantilever  Beam 


The  process  was  conducted  in  the  same  fashion  as  the  mass  spring  system.  The 
major  differences  are  the  number  of  holonomic  constraints  and  our  ability  to  choose  how 
many  modes  to  keep  when  coupling  the  system  from  substructure.  The  upcoming  tables 
will  illustrate  the  effects  of  keeping  only  the  first  few  modes  of  each  substructure  to 
calculate  for  the  dynamic  response  of  the  system  and  compare  that  to  a  system  solution 
check  that  was  generated  using  method  one  which  was  the  basic  method  of  coupling. 
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3. 


Plate  Structure 


In  this  section  Patran  and  Nastran  was  put  to  use  in  producing  natural  frequencies 
and  mode  shapes  for  two  plate  structures.  The  results  from  Nastran  where  then  extracted 
and  transferred  over  to  MATLAB  where  we  can  make  further  use  of  the  results  to  aid 
couple  the  two  plate  structures.  As  in  the  mass  spring  system  we  make  use  of  the 
interconnected  elements  of  each  substructure  to  create  the  transformation  matrix  that  is 
required  to  create  the  system.  Once  the  transformation  matrix  has  been  created  we  can 
now  assume  random  values  for  substructure  modal  damping  and  proceed  with  the 
sequence  of  events  of  method  two  to  reproduce  system  properties  such  as  mode  shape, 
natural  frequencies  and  modal  damping.  Furthennore,  by  creating  a  finite  element  model 
of  the  system  in  Patran  and  Nastran,  we  are  able  to  compare  the  values  of  co  from  method 
two,  to  the  solution  of  the  finite  element;  this  then  gave  us  further  assurance  that  the 
method  of  substructure  coupling  was  indeed  correct.  The  main  motivation  in  utilizing 
commercial  software  was  to  prevent  having  to  create  a  program  in  MATLAB  that  would 
generate  a  Finite  Element  Model  of  plate  structures  which  would  have  been  very  time 
consuming. 
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IV.  RESULTS 


Nomenclature  for  Tables 

CQ  ub  =  upper  bound  for  initial  damping  for  each  substructure 
C,a  lb  =  lower  bound  for  initial  damping  for  each  substructure 
C,subs  ub  =  upper  bound  for  substructure  damping 
C,subs  lb  =  lower  bound  for  substructure  damping 

C,subs max  =  max  value  of  substructure  damping  that  satisfies  both  substructure  and 
system  upper  bounds 

C,subs  min  =  min  value  of  substructure  damping  that  satisfies  both  substructure  and 
system  upper  bounds 

Cm  ub  =  upper  bound  for  desired  system  damping 
£  lb  =  lower  bound  for  desired  system  damping 

£  max  =  max  value  of  system  damping  that  satisfies  system  upper  bound 
Q  min  =  min  value  of  system  damping  that  satisfies  system  upper  bound 
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A 


TABLES 


Check  Wn 

Check  Zeta 

Method  3  Wn 

Method  3  Zeta 

Error  Wn  Check  Vs  Method  3 

Error  Zeta  Check  Vs  Method  3 

14.0032 

0.05 

14.00315 

0.05 

0.00% 

0.00% 

38.6003 

0.04967 

38.60027 

0.0497 

0.00% 

0.00% 

75.6722 

0.05 

75.67218 

0.05 

0.00% 

0.00% 

125.091 

0.04965 

125.0908 

0.0497 

0.00% 

0.00% 

186.867 

0.05 

186.8666 

0.05 

0.00% 

0.00% 

261.001 

0.04983 

261.0014 

0.0498 

0.00% 

0.00% 

347.5 

0.05 

347.4996 

0.05 

0.00% 

0.00% 

446.368 

0.04989 

446.3678 

0.0499 

0.00% 

0.00% 

557.616 

0.05 

557.6159 

0.05 

0.00% 

0.00% 

681.258 

0.04992 

681.2576 

0.0499 

0.00% 

0.00% 

817.312 

0.05 

817.3119 

0.05 

0.00% 

0.00% 

965.804 

0.04994 

965.8035 

0.0499 

0.00% 

0.00% 

1126.76 

0.05 

1126.764 

0.05 

0.00% 

0.00% 

1300.23 

0.04994 

1300.234 

0.0499 

0.00% 

0.00% 

1486.26 

0.05 

1486.261 

0.05 

0.00% 

0.00% 

1684.91 

0.04995 

1684.906 

0.0499 

0.00% 

0.00% 

1896.24 

0.05 

1896.24 

0.05 

0.00% 

0.00% 

2120.35 

0.04995 

2120.346 

0.0499 

0.00% 

0.00% 

2357.32 

0.05 

2357.322 

0.05 

0.00% 

0.00% 

2607.28 

0.04995 

2607.28 

0.0499 

0.00% 

0.00% 

2870.35 

0.05 

2870.348 

0.05 

0.00% 

0.00% 

3146.67 

0.04995 

3146.672 

0.0499 

0.00% 

0.00% 

3436.41 

0.05 

3436.412 

0.05 

0.00% 

0.00% 

3739.75 

0.04995 

3739.747 

0.0499 

0.00% 

0.00% 

4056.87 

0.05 

4056.873 

0.05 

0.00% 

0.00% 

4388 

0.04995 

4387.999 

0.05 

0.00% 

0.00% 

4733.35 

0.05 

4733.348 

0.05 

0.00% 

0.00% 

5093.15 

0.04996 

5093.151 

0.05 

0.00% 

0.00% 

5467.64 

0.05 

5467.636 

0.05 

0.00% 

0.00% 

5857.02 

0.04997 

5857.019 

0.05 

0.00% 

0.00% 

6261.47 

0.05 

6261.473 

0.05 

0.00% 

0.00% 

6681.09 

0.05 

6681.085 

0.05 

0.00% 

0.00% 

7115.77 

0.05 

7115.767 

0.05 

0.00% 

0.00% 

7565.08 

0.05 

7565.077 

0.05 

0.00% 

0.00% 

8027.82 

0.04994 

8027.819 

0.0499 

0.00% 

0.00% 

8501.04 

0.05 

8501.042 

0.05 

0.00% 

0.00% 

8977.19 

0.0498 

8977.192 

0.0498 

0.00% 

0.00% 

9434.6 

0.05 

9434.603 

0.05 

0.00% 

0.00% 

9806.41 

0.04955 

9806.414 

0.0496 

0.00% 

0.00% 

11155.6 

0.04907 

11155.64 

0.0491 

0.00% 

0.00% 

11598.2 

0.05 

11598.25 

0.05 

0.00% 

0.00% 

12153.2 

0.04934 

12153.21 

0.0493 

0.00% 

0.00% 

12764.6 

0.05 

12764.62 

0.05 

0.00% 

0.00% 

13415.1 

0.04946 

13415.1 

0.0495 

0.00% 

0.00% 

14098.8 

0.05 

14098.79 

0.05 

0.00% 

0.00% 

14813.7 

0.0495 

14813.73 

0.0495 

0.00% 

0.00% 

15559.5 

0.05 

15559.54 

0.05 

0.00% 

0.00% 

16336.6 

0.04951 

16336.58 

0.0495 

0.00% 

0.00% 

17145.6 

0.05 

17145.62 

0.05 

0.00% 

0.00% 

17987.7 

0.04949 

17987.65 

0.0495 

0.00% 

0.00% 

18863.8 

0.05 

18863.84 

0.05 

0.00% 

0.00% 

19775.4 

0.04945 

19775.41 

0.0495 

0.00% 

0.00% 

20723.7 

0.05 

20723.67 

0.05 

0.00% 

0.00% 

21709.9 

0.0494 

21709.92 

0.0494 

0.00% 

0.00% 

22735.4 

0.05 

22735.44 

0.05 

0.00% 

0.00% 
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23801.4 

0.04932 

23801.43 

0.0493 

0.00% 

0.00% 

24908.9 

0.05 

24908.95 

0.05 

0.00% 

0.00% 

26058.8 

0.04922 

26058.84 

0.0492 

0.00% 

0.00% 

27251.6 

0.05 

27251.65 

0.05 

0.00% 

0.00% 

28487.5 

0.04911 

28487.47 

0.0491 

0.00% 

0.00% 

29765.9 

0.05 

29765.85 

0.05 

0.00% 

0.00% 

31085.6 

0.04898 

31085.57 

0.049 

0.00% 

0.00% 

32444.4 

0.05 

32444.44 

0.05 

0.00% 

0.00% 

33839.1 

0.04885 

33839.07 

0.0489 

0.00% 

0.00% 

35264.6 

0.05 

35264.56 

0.05 

0.00% 

0.00% 

36714.2 

0.04879 

36714.16 

0.0488 

0.00% 

0.00% 

38179 

0.05 

38179 

0.05 

0.00% 

0.00% 

39647.7 

0.04893 

39647.66 

0.0489 

0.00% 

0.00% 

41106 

0.05 

41105.96 

0.05 

0.00% 

0.00% 

42536.7 

0.04955 

42536.7 

0.0496 

0.00% 

0.00% 

43919.6 

0.05 

43919.6 

0.05 

0.00% 

0.00% 

45231.5 

0.05 

45231.53 

0.05 

0.00% 

0.00% 

46447 

0.04863 

46446.97 

0.0486 

0.00% 

0.00% 

47538.9 

0.05 

47538.86 

0.05 

0.00% 

0.00% 

48479.8 

0.04443 

48479.83 

0.0444 

0.00% 

0.00% 

49243.8 

0.04931 

49243.85 

0.0493 

0.00% 

0.00% 

49808 

0.0347 

49807.96 

0.0347 

0.00% 

0.00% 

50154.1 

0.05 

50154.14 

0.05 

0.00% 

0.00% 

Table  1.  20  Element  Per  Substructure  All  modes  Kept 
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Check  Wn 

Check  Zeta 

Method  3  Wn 

Method  3  Zeta 

Error  Wn  Check  Vs  Method  3 

Error  Zeta  Check  Vs  Method  3 

14.0032 

0.05 

14.09988 

0.05 

0.69% 

0.00% 

38.6003 

0.04967 

38.6004 

0.0497 

0.00% 

0.12% 

75.6722 

0.05 

76.3976 

0.05 

0.96% 

0.00% 

125.091 

0.04965 

125.0923 

0.0497 

0.00% 

0.06% 

186.867 

0.05 

188.6633 

0.05 

0.96% 

0.00% 

261.001 

0.04983 

261.0077 

0.0498 

0.00% 

0.03% 

347.5 

0.05 

350.8928 

0.05 

0.98% 

0.00% 

446.368 

0.04989 

446.3863 

0.0499 

0.00% 

0.02% 

557.616 

0.05 

563.1438 

0.05 

0.99% 

0.00% 

681.258 

0.04992 

681.3006 

0.0499 

0.01% 

0.01% 

817.312 

0.05 

825.5409 

0.05 

1.01% 

0.00% 

965.804 

0.04994 

965.89 

0.0499 

0.01% 

0.01% 

1126.76 

0.05 

1138.291 

0.05 

1 .02% 

0.00% 

1300.23 

0.04994 

1300.39 

0.05 

0.01% 

0.01% 

1486.26 

0.05 

1501.721 

0.05 

1 .04% 

0.00% 

1684.91 

0.04995 

1685.169 

0.05 

0.02% 

0.02% 

1896.24 

0.05 

1916.312 

0.05 

1 .06% 

0.00% 

2120.35 

0.04995 

2120.763 

0.05 

0.02% 

0.03% 

2357.32 

0.05 

2382.742 

0.05 

1 .08% 

0.00% 

2607.28 

0.04995 

2607.91 

0.05 

0.02% 

0.04% 

2870.35 

0.05 

2901.926 

0.05 

1.10% 

0.00% 

3146.67 

0.04995 

3147.59 

0.05 

0.03% 

0.06% 

3436.41 

0.05 

3475.054 

0.05 

1.12% 

0.00% 

3739.75 

0.04995 

3741 .042 

0.05 

0.03% 

0.08% 

4056.87 

0.05 

4103.614 

0.05 

1.15% 

0.00% 

4388 

0.04995 

4389.776 

0.05 

0.04% 

0.10% 

4733.35 

0.05 

4789.407 

0.05 

1.18% 

0.01% 

5093.15 

0.04996 

5095.533 

0.05 

0.05% 

0.09% 

5467.64 

0.05 

5534.501 

0.05 

1 .22% 

0.04% 

5857.02 

0.04997 

5860.137 

0.05 

0.05% 

i  MEP 

6261.47 

0.05 

6341.06 

0.05 

1.27% 

0.10% 

6681.09 

0.05 

6685.061 

0.05 

0.06% 

0.00% 

7115.77 

0.05 

7210.761 

0.0499 

1 .33% 

0.18% 

7565.08 

0.05 

7569.964 

0.05 

0.06% 

0.00% 

8027.82 

0.04994 

8142.548 

0.0498 

1 .43% 

0.20% 

8501.04 

0.05 

8506.581 

0.05 

0.07% 

0.00% 

8977.19 

0.0498 

9120.192 

0.0497 

1 .59% 

0.24% 

9434.6 

0.05 

9439.023 

0.05 

0.05% 

0.00% 

9806.41 

0.04955 

10043.56 

0.0489 

2.42% 

1.41% 

11155.6 

0.04907 

11174.97 

0.0493 

0.52% 

11598.2 

0.05 

11685.07 

0.05 

12153.2 

0.04934 

12173.57 

0.0496 

0.17% 

0.47% 

12764.6 

0.05 

12906.03 

0.05 

1.11% 

0.00% 

13415.1 

0.04946 

13438.44 

0.0497 

0.17% 

0.50% 

14098.8 

0.05 

14279.65 

0.05 

1 .28% 

0.00% 

14813.7 

0.0495 

14841.67 

0.0498 

0.19% 

0.56% 

15559.5 

0.05 

15779.88 

0.05 

1 .42% 

0.00% 

16336.6 

0.04951 

16370.63 

0.0498 

0.21% 

0.63% 

17145.6 

0.05 

17410.08 

0.05 

1 .54% 

0.00% 

17987.7 

0.04949 

18029.6 

0.0498 

0.23% 

0.72% 

18863.8 

0.05 

19179.88 

0.05 

1 .68% 

0.00% 

19775.4 

0.04945 

19827.49 

0.0499 

0.26% 

0.83% 

20723.7 

0.05 

21101.34 

0.05 

0.00% 

21709.9 

0.0494 

21775.06 

0.0499 

22735.4 

0.05 

23187.71 

0.05 

23801.4 

0.04932 

23883.55 

0.0499 

24908.9 

0.05 

25452.22 

0.05 

26058.8 

0.04922 

26163.25 

0.0499 

0.40% 

1.37% 

27251.6 

0.05 

27906.36 

0.05 

2.40% 

0.00% 

28487.5 

0.04911 

28621.42 

0.0499 

0.47% 

1.61% 

29765.9 

0.05 

30557.1 

0.05 

2.66% 

0.00% 

31085.6 

0.04898 

31258.96 

0.0499 

0.56% 

1 .90% 

32444.4 

0.05 

33402.51 

0.05 

2.95% 

0.00% 

33839.1 

0.04885 

34065.49 

0.0499 

r ;vr.' 

2.16% 

35264.6 

0.05 

36426.07 

0.05 

0.00% 

36714.2 

0.04879 

37012.26 

0.0499 

0.81% 

2.30% 

38179 

0.05 

39596.4 

0.05 

3.71% 

0.00% 

39647.7 

0.04893 

40045.74 

0.0499 

41106 

0.05 

42536.7 

0.04955 

43919.6 

0.05 

45231.5 

0.05 

46447 

0.04863 

47538.9 

0.05 

48479.8 

0.04443 

49243.8 

0.04931 

49808 

0.0347 

50154.1 

0.05 

Table  2.  20  Element  Per  Substructure  35  modes  Kept 
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Check  Wn 

Check  Zeta 

Method  3  Wn 

Method  3  Zeta 

Error  Wn  Check  Vs  Method  3 

Error  Zeta  Check  Vs  Method  3 

14.0032 

0.05 

14.12993 

0.05 

0.91% 

0.00% 

38.6003 

0.04967 

38.60047 

0.0497 

0.00% 

0.16% 

75.6722 

0.05 

76.62695 

0.05 

1.26% 

0.00% 

125.091 

0.04965 

125.0931 

0.0497 

0.00% 

0.07% 

186.867 

0.05 

189.2424 

0.05 

1.27% 

0.00% 

261.001 

0.04983 

261.0113 

0.0498 

0.00% 

0.02% 

347.5 

0.05 

352.0073 

0.05 

1.30% 

0.00% 

446.368 

0.04989 

446.3968 

0.0499 

0.01% 

0.00% 

557.616 

0.05 

564.9938 

0.05 

1.32% 

0.00% 

681.258 

0.04992 

681.3252 

0.0499 

0.01% 

0.02% 

817.312 

0.05 

828.3466 

0.05 

1 .35% 

0.00% 

965.804 

0.04994 

965.9393 

0.0499 

0.01% 

0.03% 

1126.76 

0.05 

1142.295 

0.05 

1 .38% 

0.00% 

1300.23 

0.04994 

1300.48 

0.0499 

0.02% 

0.03% 

1486.26 

0.05 

1507.19 

0.05 

1.41% 

0.00% 

1684.91 

0.04995 

1685.319 

0.0499 

0.02% 

0.03% 

1896.24 

0.05 

1923.543 

0.05 

1.44% 

0.00% 

2120.35 

0.04995 

2121 

0.0499 

0.03% 

0.03% 

2357.32 

0.05 

2392.067 

0.05 

1.47% 

0.00% 

2607.28 

0.04995 

2608.269 

0.0499 

0.04% 

0.02% 

2870.35 

0.05 

2913.721 

0.05 

1.51% 

0.00% 

3146.67 

0.04995 

3148.112 

0.0499 

0.05% 

0.01% 

3436.41 

0.05 

3489.746 

0.05 

1 .55% 

0.00% 

3739.75 

0.04995 

3741.78 

0.0499 

0.05% 

0.00% 

4056.87 

0.05 

4121.701 

0.05 

1.60% 

0.00% 

4388 

0.04995 

4390.79 

0.05 

0.06% 

0.01% 

4733.35 

0.05 

4811.474 

0.05 

1 .65% 

0.00% 

5093.15 

0.04996 

5096.892 

0.05 

0.07% 

0.02% 

5467.64 

0.05 

5561.256 

0.05 

1.71% 

0.00% 

5857.02 

0.04997 

5861.919 

0.05 

0.08% 

0.06% 

6261.47 

0.05 

6373.383 

0.05 

1.79% 

0.00% 

6681.09 

0.05 

6687.34 

0.05 

0.09% 

0.00% 

7115.77 

0.05 

7249.799 

0.05 

1 .88% 

0.06% 

7565.08 

0.05 

7572.772 

0.05 

0.10% 

0.00% 

8027.82 

0.04994 

8189.856 

0.0499 

2.02% 

0.06% 

8501.04 

0.05 

8509.775 

0.05 

0.10% 

0.00% 

8977.19 

0.0498 

9177.231 

0.0498 

2.23% 

0.02% 

9434.6 

0.05 

9441.573 

0.05 

0.07% 

0.00% 

9806.41 

0.04955 

10190.56 

0.0488 

3.92% 

1.54% 

11155.6 

0.04907 

11185.92 

0.0494 

0.27% 

0.70% 

11598.2 

0.05 

11732.62 

0.05 

1.16% 

0.00% 

12153.2 

0.04934 

12185.58 

0.0496 

0.27% 

0.59% 

12764.6 

0.05 

12979.25 

0.05 

1 .68% 

0.00% 

13415.1 

0.04946 

13452.45 

0.0497 

0.28% 

0.42% 

14098.8 

0.05 

14373.69 

0.05 

1 .95% 

0.00% 

14813.7 

0.0495 

14858.75 

0.0499 

0.30% 

0.79% 

15559.5 

0.05 

15897.01 

0.05 

2.17% 

0.00% 

16336.6 

0.04951 

16391.98 

0.0497 

0.34% 

0.31% 

17145.6 

0.05 

17555.52 

0.05 

2.39% 

0.00% 

17987.7 

0.04949 

18056.78 

0.05 

0.38% 

1 .03% 

18863.8 

0.05 

19361.93 

0.05 

2.64% 

0.00% 

19775.4 

0.04945 

19862.83 

0.0497 

0.44% 

0.58% 

20723.7 

0.05 

21333.36 

0.05 

2.94% 

0.00% 

21709.9 

0.0494 

21822.38 

0.0499 

0.52% 

1.07% 

22735.4 

0.05 

23493.38 

0.05 

3.33% 

0.00% 

23801.4 

0.04932 

23950.06 

0.0499 

0.62% 

1.10% 

24908.9 

0.05 

25885.77 

0.05 

3.92% 

0.00% 

26058.8 

0.04922 

26267.28 

0.0499 

0.80% 

1 .36% 

27251.6 

0.05 

28487.5 

0.04911 

29765.9 

0.05 

31085.6 

0.04898 

32444.4 

0.05 

33839.1 

0.04885 

35264.6 

0.05 

36714.2 

0.04879 

38179 

0.05 

39647.7 

0.04893 

41106 

0.05 

42536.7 

0.04955 

43919.6 

0.05 

45231.5 

0.05 

46447 

0.04863 

47538.9 

0.05 

48479.8 

0.04443 

49243.8 

0.04931 

49808 

0.0347 

50154.1 

0.05 

Table  3.  20  Element  Per  Substructure  30  modes  Kept 
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Check  Wn 

Check  Zeta 

Method  3  Wn 

Method  3  Zeta 

Error  Wn  Check  Vs  Method  3 

Error  Zeta  Check  Vs  Method  3 

14.0032 

0.05 

14.17136 

0.05 

1.20% 

0.00% 

38.6003 

0.04967 

38.60066 

0.0498 

0.00% 

0.20% 

75.6722 

0.05 

76.9464 

0.05 

1 .68% 

0.00% 

125.091 

0.04965 

125.0952 

0.0497 

0.00% 

0.07% 

186.867 

0.05 

190.0578 

0.05 

1.71% 

0.00% 

261.001 

0.04983 

261.0204 

0.0498 

0.01% 

0.00% 

347.5 

0.05 

353.5932 

0.05 

1.75% 

0.00% 

446.368 

0.04989 

446.4232 

0.0499 

0.01% 

0.04% 

557.616 

0.05 

567.6543 

0.05 

1 .80% 

0.00% 

681.258 

0.04992 

681.3866 

0.0499 

0.02% 

0.08% 

817.312 

0.05 

832.4239 

0.05 

1 .85% 

0.00% 

965.804 

0.04994 

966.0627 

0.0499 

0.03% 

0.10% 

1126.76 

0.05 

1148.174 

0.05 

1 .90% 

0.00% 

1300.23 

0.04994 

1300.703 

0.0499 

0.04% 

0.11% 

1486.26 

0.05 

1515.305 

0.05 

1 .95% 

0.00% 

1684.91 

0.04995 

1685.695 

0.0499 

0.05% 

0.12% 

1896.24 

0.05 

1934.384 

0.05 

2.01% 

0.00% 

2120.35 

0.04995 

2121.595 

0.0499 

0.06% 

0.13% 

2357.32 

0.05 

2406.196 

0.05 

2.07% 

0.00% 

2607.28 

0.04995 

2609.17 

0.0499 

0.07% 

0.13% 

2870.35 

0.05 

2931.781 

0.05 

2.14% 

0.00% 

3146.67 

0.04995 

3149.426 

0.0499 

0.09% 

0.13% 

3436.41 

0.05 

3512.486 

0.05 

2.21% 

0.00% 

3739.75 

0.04995 

3743.639 

0.0499 

0.10% 

0.12% 

4056.87 

0.05 

4149.998 

0.05 

2.30% 

0.00% 

4388 

0.04995 

4393.354 

0.0499 

0.12% 

0.12% 

4733.35 

0.05 

4846.373 

0.05 

2.39% 

0.00% 

5093.15 

0.04996 

5100.35 

0.0499 

0.14% 

0.11% 

5467.64 

0.05 

5604.01 

0.05 

2.49% 

0.00% 

5857.02 

0.04997 

5866.485 

0.0499 

0.16% 

0.12% 

6261.47 

0.05 

6425.508 

0.05 

2.62% 

0.00% 

6681.09 

0.05 

6693.23 

0.0499 

0.18% 

0.14% 

7115.77 

0.05 

7313.077 

0.05 

2.77% 

0.00% 

7565.08 

0.05 

7580.122 

0.0499 

0.20% 

0.11% 

8027.82 

0.04994 

8265.936 

0.05 

2.97% 

0.12% 

8501.04 

0.05 

8518.26 

0.05 

0.20% 

0.00% 

8977.19 

0.0498 

9262.239 

0.05 

3.18% 

0.31% 

9434.6 

0.05 

9448.436 

0.05 

0.15% 

0.00% 

9806.41 

0.04955 

10510.54 

0.049 

7.18% 

1 .02% 

11155.6 

0.04907 

11215.53 

0.0498 

0.54% 

1 .47% 

11598.2 

0.05 

11846.5 

0.05 

2.14% 

0.00% 

12153.2 

0.04934 

12220.83 

0.0495 

0.56% 

0.24% 

12764.6 

0.05 

13147.17 

0.05 

3.00% 

0.00% 

13415.1 

0.04946 

13496.48 

0.05 

0.61% 

1 .09% 

14098.8 

0.05 

14601.45 

0.05 

3.57% 

0.00% 

14813.7 

0.0495 

14918.37 

0.0497 

0.71% 

0.49% 

15559.5 

0.05 

16222.34 

0.05 

4.26% 

0.00% 

16336.6 

0.04951 

16482.55 

0.0499 

0.89% 

0.75% 

17145.6 

0.05 

17987.7 

0.04949 

18863.8 

0.05 

19775.4 

0.04945 

20723.7 

0.05 

21709.9 

0.0494 

22735.4 

0.05 

23801.4 

0.04932 

24908.9 

0.05 

26058.8 

0.04922 

27251.6 

0.05 

28487.5 

0.04911 

29765.9 

0.05 

31085.6 

0.04898 

32444.4 

0.05 

33839.1 

0.04885 

35264.6 

0.05 

36714.2 

0.04879 

38179 

0.05 

39647.7 

0.04893 

41106 

0.05 

42536.7 

0.04955 

43919.6 

0.05 

45231.5 

0.05 

46447 

0.04863 

47538.9 

0.05 

48479.8 

0.04443 

49243.8 

0.04931 

49808 

0.0347 

50154.1 

0.05 

Table  4.  20  Element  Per  Substructure  25  modes  Kept 
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Check  Wn 

Check  Zeta 

Method  3  Wn 

Method  3  Zeta 

Error  Wn  Check  Vs  Method  3 

Error  Zeta  Check  Vs  Method  3 

14.0032 

0.05 

14.21897 

0.05 

1 .54% 

0.00% 

38.6003 

0.04967 

38.60105 

0.0498 

0.00% 

0.24% 

75.6722 

0.05 

77.31798 

0.05 

2.17% 

0.00% 

125.091 

0.04965 

125.0995 

0.0497 

0.01% 

0.06% 

186.867 

0.05 

191.0191 

0.05 

2.22% 

0.00% 

261.001 

0.04983 

261.0391 

0.0498 

0.01% 

0.04% 

347.5 

0.05 

355.4874 

0.05 

2.30% 

0.00% 

446.368 

0.04989 

446.4781 

0.0498 

0.02% 

0.11% 

557.616 

0.05 

570.8732 

0.05 

2.38% 

0.00% 

681.258 

0.04992 

681.5146 

0.0498 

0.04% 

0.16% 

817.312 

0.05 

837.4212 

0.05 

2.46% 

0.00% 

965.804 

0.04994 

966.32 

0.0498 

0.05% 

0.19% 

1126.76 

0.05 

1155.473 

0.05 

2.55% 

0.00% 

1300.23 

0.04994 

1301.17 

0.0498 

0.07% 

0.22% 

1486.26 

0.05 

1525.511 

0.05 

2.64% 

0.00% 

1684.91 

0.04995 

1686.48 

0.0498 

0.09% 

0.24% 

1896.24 

0.05 

1948.202 

0.05 

2.74% 

0.00% 

2120.35 

0.04995 

2122.842 

0.0498 

0.12% 

0.25% 

2357.32 

0.05 

2424.451 

0.05 

2.85% 

0.00% 

2607.28 

0.04995 

2611.063 

0.0498 

0.15% 

0.26% 

2870.35 

0.05 

2955.449 

0.05 

2.96% 

0.00% 

3146.67 

0.04995 

3152.203 

0.0498 

0.18% 

0.26% 

3436.41 

0.05 

3542.739 

0.05 

3.09% 

0.00% 

3739.75 

0.04995 

3747.599 

0.0498 

0.21% 

0.26% 

4056.87 

0.05 

4188.262 

0.05 

3.24% 

0.00% 

4388 

0.04995 

4398.88 

0.0498 

0.25% 

0.25% 

4733.35 

0.05 

4894.407 

0.05 

3.40% 

0.00% 

5093.15 

0.04996 

5107.925 

0.0498 

0.29% 

0.27% 

5467.64 

0.05 

5664.005 

0.05 

3.59% 

0.00% 

5857.02 

0.04997 

5876.723 

0.0498 

0.34% 

0.27% 

6261.47 

0.05 

6500.173 

0.05 

3.81% 

0.00% 

6681.09 

0.05 

6706.887 

0.0498 

0.39% 

0.32% 

7115.77 

0.05 

7405.475 

0.05 

4.07% 

0.00% 

7565.08 

0.05 

7597.999 

0.0499 

0.44% 

0.29% 

8027.82 

0.04994 

8377.479 

0.05 

4.36% 

0.12% 

8501.04 

0.05 

8540.406 

0.0499 

0.46% 

0.24% 

8977.19 

0.0498 

9373.524 

0.05 

4.41% 

0.40% 

9434.6 

0.05 

9468.096 

0.0499 

0.36% 

0.14% 

9806.41 

0.04955 

11155.6 

0.04907 

11598.2 

0.05 

12153.2 

0.04934 

12764.6 

0.05 

13415.1 

0.04946 

14098.8 

0.05 

14813.7 

0.0495 

15559.5 

0.05 

16336.6 

0.04951 

17145.6 

0.05 

17987.7 

0.04949 

18863.8 

0.05 

19775.4 

0.04945 

20723.7 

0.05 

21709.9 

0.0494 

22735.4 

0.05 

23801.4 

0.04932 

24908.9 

0.05 

26058.8 

0.04922 

27251.6 

0.05 

28487.5 

0.04911 

29765.9 

0.05 

31085.6 

0.04898 

32444.4 

0.05 

33839.1 

0.04885 

35264.6 

0.05 

36714.2 

0.04879 

38179 

0.05 

39647.7 

0.04893 

41106 

0.05 

42536.7 

0.04955 

43919.6 

0.05 

45231.5 

0.05 

46447 

0.04863 

47538.9 

0.05 

48479.8 

0.04443 

49243.8 

0.04931 

49808 

0.0347 

50154.1 

0.05 

Table  5.  20  Element  Per  Substructure  20  modes  Kept 
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Check  Wn 

Check  Zeta 

Method  3  Wn 

Method  3  Zeta 

Error  Wn  Check  Vs  Method  3 

Error  Zeta  Check  Vs  Method  3 

14.0032 

0.05 

14.2942 

0.05 

2.08% 

0.00% 

38.6003 

0.04967 

38.60231 

0.0498 

0.01% 

0.29% 

75.6722 

0.05 

77.91517 

0.05 

2.96% 

0.00% 

125.091 

0.04965 

125.1135 

0.0497 

0.02% 

0.04% 

186.867 

0.05 

192.5928 

0.05 

3.06% 

0.00% 

261.001 

0.04983 

261.0998 

0.0498 

0.04% 

0.11% 

347.5 

0.05 

358.6457 

0.05 

3.21% 

0.00% 

446.368 

0.04989 

446.6556 

0.0498 

0.06% 

0.20% 

557.616 

0.05 

576.3424 

0.05 

3.36% 

0.00% 

681.258 

0.04992 

681.9285 

0.0498 

0.10% 

0.27% 

817.312 

0.05 

846.0811 

0.05 

3.52% 

0.00% 

965.804 

0.04994 

967.1543 

0.0498 

0.14% 

0.31% 

1126.76 

0.05 

1168.397 

0.05 

3.69% 

0.00% 

1300.23 

0.04994 

1302.692 

0.0498 

0.19% 

0.34% 

1486.26 

0.05 

1544.029 

0.05 

3.89% 

0.00% 

1684.91 

0.04995 

1689.065 

0.0498 

0.25% 

0.35% 

1896.24 

0.05 

1974.008 

0.05 

4.10% 

0.00% 

2120.35 

0.04995 

2127.016 

0.0498 

0.31% 

0.36% 

2357.32 

0.05 

2459.782 

0.05 

4.35% 

0.00% 

2607.28 

0.04995 

2617.57 

0.0498 

0.39% 

0.35% 

2870.35 

0.05 

3003.435 

0.05 

4.64% 

0.00% 

3146.67 

0.04995 

3162.139 

0.0498 

0.49% 

0.33% 

3436.41 

0.05 

3608.124 

0.05 

5.00% 

0.00% 

3739.75 

0.04995 

3762.706 

0.0498 

0.61% 

0.31% 

4056.87 

0.05 

4279.266 

0.05 

5.48% 

0.00% 

4388 

0.04995 

4422.318 

0.0498 

0.78% 

0.29% 

4733.35 

0.05 

5029.64 

0.05 

6.26% 

0.00% 

5093.15 

0.04996 

5147.391 

0.0498 

1 .06% 

0.25% 

5467.64 

0.05 

5857.02 

0.04997 

6261.47 

0.05 

6681.09 

0.05 

7115.77 

0.05 

7565.08 

0.05 

8027.82 

0.04994 

8501.04 

0.05 

8977.19 

0.0498 

9434.6 

0.05 

9806.41 

0.04955 

11155.6 

0.04907 

11598.2 

0.05 

12153.2 

0.04934 

12764.6 

0.05 

13415.1 

0.04946 

14098.8 

0.05 

14813.7 

0.0495 

15559.5 

0.05 

16336.6 

0.04951 

17145.6 

0.05 

17987.7 

0.04949 

18863.8 

0.05 

19775.4 

0.04945 

20723.7 

0.05 

21709.9 

0.0494 

22735.4 

0.05 

23801.4 

0.04932 

24908.9 

0.05 

26058.8 

0.04922 

27251.6 

0.05 

28487.5 

0.04911 

29765.9 

0.05 

31085.6 

0.04898 

32444.4 

0.05 

33839.1 

0.04885 

35264.6 

0.05 

36714.2 

0.04879 

38179 

0.05 

39647.7 

0.04893 

41106 

0.05 

42536.7 

0.04955 

43919.6 

0.05 

45231.5 

0.05 

46447 

0.04863 

47538.9 

0.05 

48479.8 

0.04443 

49243.8 

0.04931 

49808 

0.0347 

50154.1 

0.05 

Table  6.  20  Element  Per  Substructure  15  modes  Kept 
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Check  Wn 

Check  Zeta 

Method  3  Wn 

Method  3  Zeta 

Error  Wn  Check  Vs  Method  3 

Error  Zeta  Check  Vs  Method  3 

14.0032 

0.05 

14.45803 

0.05 

3.25% 

0.00% 

38.6003 

0.04967 

38.60776 

0.0498 

0.02% 

0.36% 

75.6722 

0.05 

79.25888 

0.05 

4.74% 

0.00% 

125.091 

0.04965 

125.1741 

0.0496 

0.07% 

0.07% 

186.867 

0.05 

196.2635 

0.05 

5.03% 

0.00% 

261.001 

0.04983 

261.3633 

0.0497 

0.14% 

0.28% 

347.5 

0.05 

366.2861 

0.05 

5.41% 

0.00% 

446.368 

0.04989 

447.4288 

0.0496 

0.24% 

0.49% 

557.616 

0.05 

590.1081 

0.05 

5.83% 

0.00% 

681.258 

0.04992 

683.7469 

0.0497 

0.37% 

0.35% 

817.312 

0.05 

868.9125 

0.05 

6.31% 

0.00% 

965.804 

0.04994 

970.8904 

0.0495 

0.53% 

0.78% 

1126.76 

0.05 

1204.563 

0.05 

6.90% 

0.00% 

1300.23 

0.04994 

1309.773 

0.0498 

0.73% 

0.24% 

1486.26 

0.05 

1600.535 

0.05 

7.69% 

0.00% 

1684.91 

0.04995 

1702.045 

0.0496 

1 .02% 

0.64% 

1896.24 

0.05 

2065.763 

0.05 

8.94% 

0.00% 

2120.35 

0.04995 

2151.829 

0.0498 

1 .48% 

0.36% 

2357.32 

0.05 

2607.28 

0.04995 

2870.35 

0.05 

3146.67 

0.04995 

3436.41 

0.05 

3739.75 

0.04995 

4056.87 

0.05 

4388 

0.04995 

4733.35 

0.05 

5093.15 

0.04996 

5467.64 

0.05 

5857.02 

0.04997 

6261.47 

0.05 

6681.09 

0.05 

7115.77 

0.05 

7565.08 

0.05 

8027.82 

0.04994 

8501.04 

0.05 

8977.19 

0.0498 

9434.6 

0.05 

9806.41 

0.04955 

11155.6 

0.04907 

11598.2 

0.05 

12153.2 

0.04934 

12764.6 

0.05 

13415.1 

0.04946 

14098.8 

0.05 

14813.7 

0.0495 

15559.5 

0.05 

16336.6 

0.04951 

17145.6 

0.05 

17987.7 

0.04949 

18863.8 

0.05 

19775.4 

0.04945 

20723.7 

0.05 

21709.9 

0.0494 

22735.4 

0.05 

23801.4 

0.04932 

24908.9 

0.05 

26058.8 

0.04922 

27251.6 

0.05 

28487.5 

0.04911 

29765.9 

0.05 

31085.6 

0.04898 

32444.4 

0.05 

33839.1 

0.04885 

35264.6 

0.05 

36714.2 

0.04879 

38179 

0.05 

39647.7 

0.04893 

41106 

0.05 

42536.7 

0.04955 

43919.6 

0.05 

45231.5 

0.05 

46447 

0.04863 

47538.9 

0.05 

48479.8 

0.04443 

49243.8 

0.04931 

49808 

0.0347 

50154.1 

0.05 

Table  7.  20  Element  Per  Substructure  10  modes  Kept 
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Check  Wn 

Check  Zeta 

Method  3  Wn 

Method  3  Zeta 

Error  Wn  Check  Vs  Method  3 

Error  Zeta  Check  Vs  Method  3 

14.0032 

0.05 

15.00189 

0.05 

7.13% 

0.00% 

38.6003 

0.04967 

38.66069 

0.0497 

0.16% 

0.02% 

75.6722 

0.05 

84.17335 

0.05 

11.23% 

0.00% 

125.091 

0.04965 

125.7672 

0.05 

0.54% 

0.70% 

186.867 

0.05 

211.4017 

0.05 

13.13% 

0.00% 

261.001 

0.04983 

264.0277 

0.0485 

1.16% 

2.62% 

347.5 

0.05 

403.839 

0.05 

16.21% 

0.00% 

446.368 

0.04989 

456.284 

0.05 

2.22% 

0.19% 

557.616 

0.05 

681.258 

0.04992 

817.312 

0.05 

965.804 

0.04994 

1126.76 

0.05 

1300.23 

0.04994 

1486.26 

0.05 

1684.91 

0.04995 

1896.24 

0.05 

2120.35 

0.04995 

2357.32 

0.05 

2607.28 

0.04995 

2870.35 

0.05 

3146.67 

0.04995 

3436.41 

0.05 

3739.75 

0.04995 

4056.87 

0.05 

4388 

0.04995 

4733.35 

0.05 

5093.15 

0.04996 

5467.64 

0.05 

5857.02 

0.04997 

6261.47 

0.05 

6681.09 

0.05 

7115.77 

0.05 

7565.08 

0.05 

8027.82 

0.04994 

8501.04 

0.05 

8977.19 

0.0498 

9434.6 

0.05 

9806.41 

0.04955 

11155.6 

0.04907 

11598.2 

0.05 

12153.2 

0.04934 

12764.6 

0.05 

13415.1 

0.04946 

14098.8 

0.05 

14813.7 

0.0495 

15559.5 

0.05 

16336.6 

0.04951 

17145.6 

0.05 

17987.7 

0.04949 

18863.8 

0.05 

19775.4 

0.04945 

20723.7 

0.05 

21709.9 

0.0494 

22735.4 

0.05 

23801.4 

0.04932 

24908.9 

0.05 

26058.8 

0.04922 

27251.6 

0.05 

28487.5 

0.04911 

29765.9 

0.05 

31085.6 

0.04898 

32444.4 

0.05 

33839.1 

0.04885 

35264.6 

0.05 

36714.2 

0.04879 

38179 

0.05 

39647.7 

0.04893 

41106 

0.05 

42536.7 

0.04955 

43919.6 

0.05 

45231.5 

0.05 

46447 

0.04863 

47538.9 

0.05 

48479.8 

0.04443 

49243.8 

0.04931 

49808 

0.0347 

50154.1 

0.05 

Table  8.  20  Element  Per  Substructure  5  modes  Kept 
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Modes 

Kept 

C0 

^3 subs 

ub 

^3 subs 

lb 

^sys 

ub 

^>sys 

lb 

^3 subs 

max 

^3 subs 

min 

^ sys 

max 

^>sys 

min 

ALL 

0.2 

0.1 

0.01 

0.05 

0.01 

0.1 

0.0339 

0.05 

0.034705 

35 

0.2 

0.1 

0.01 

0.05 

0.01 

0.0721 

0.0501 

0.05 

0.0489 

30 

0.2 

0.1 

0.01 

0.05 

0.01 

0.0722 

0.0502 

0.05 

0.0488 

25 

0.2 

0.1 

0.01 

0.05 

0.01 

0.0723 

0.05 

0.05 

0.049 

20 

0.2 

0.1 

0.01 

0.05 

0.01 

0.0724 

0.0501 

0.05 

0.0497 

15 

0.2 

0.1 

0.01 

0.05 

0.01 

0.0725 

0.0503 

0.05 

0.0497 

10 

0.2 

0.1 

0.01 

0.05 

0.01 

0.0728 

0.0506 

0.05 

0.0495 

5 

0.2 

0.1 

0.01 

0.05 

0.01 

0.0737 

0.0521 

0.05 

0.0485 

Table  9.  For  a  20  Element  Beam  Structure  with  varying  number  of  modes  retained; 

Calculated  Results  for  Substructure  Damping  which  meets  the  prescribed  system 
damping  constraint  and  bounds  for  the  values  of  substructure  damping.  The 
results  for  system  damping  are  given  as  well. 


C0 

^3 subs 

ub 

^3 subs 

lb 

^sys 

ub 

^ sys 

lb 

max 

^3 subs 

min 

^ sys 

max 

jsys 

min 

Iterations 

Function 

Call 

0.001 

0.1 

0.01 

0.05 

0.01 

0.0724 

0.0501 

0.05 

0.0497 

2 

63 

0.005 

0.1 

0.01 

0.05 

0.01 

0.0724 

0.0501 

0.05 

0.0497 

2 

63 

0.01 

0.1 

0.01 

0.05 

0.01 

0.0724 

0.0501 

0.05 

0.0497 

2 

63 

0.1 

0.1 

0.01 

0.05 

0.01 

0.0724 

0.0501 

0.05 

0.0497 

2 

63 

0.2 

0.1 

0.01 

0.05 

0.01 

0.0724 

0.0501 

0.05 

0.0497 

2 

63 

0.5 

0.1 

0.01 

0.05 

0.01 

0.0724 

0.0501 

0.05 

0.0497 

2 

63 

1 

0.1 

0.01 

0.05 

0.01 

0.0724 

0.0501 

0.05 

0.0497 

2 

63 

10 

0.1 

0.01 

0.05 

0.01 

0.0724 

0.0501 

0.05 

0.0497 

2 

63 

Table  10.  Using  a  20  Element  Beam  Substructures  and  Keeping  only  20  modes;  System  and 
Substructure  damping  results  with  varying  start  points  (initial  values  of  the  design 

variable,  substructure  damping) 


c 

^3 subs 

ub 

^3 subs 

lb 

^sys 

ub 

^ sys 

lb 

max 

^3 subs 

min 

^ sys 

max 

jsys 

min 

Iterations 

Function 

Call 

0.2 

0.2 

0.01 

0.05 

0.01 

0.0724 

0.0501 

0.05 

0.0497 

2 

63 

0.2 

0.5 

0.01 

0.05 

0.01 

0.0724 

0.0501 

0.05 

0.0497 

2 

63 

0.2 

0.1 

0.01 

0.05 

0.01 

0.0724 

0.0501 

0.05 

0.0497 

2 

63 

0.2 

0.05 

0.01 

0.05 

0.01 

0.0500 

0.0500 

0.0498 

0.0353 

1 

42 

0.2 

0.03 

0.01 

0.05 

0.01 

0.0300 

0.0300 

0.0299 

0.0212 

1 

42 

0.2 

0.05 

0.03 

0.05 

0.01 

0.0500 

0.0500 

0.0498 

0.0353 

1 

42 

0.2 

0.04 

0.03 

0.05 

0.01 

0.0400 

0.0400 

0.0399 

0.0283 

1 

42 

0.2 

0.015 

0 

0.05 

0.01 

0.0150 

0.0150 

0.0150 

0.0106 

1 

42 

Table  11.  Using  a  20  Element  Beam  Substructures  and  Keeping  only  20  modes  and  Initial 
value  for  the  design  variable  kept  at  0.2;  System  and  Substructure  damping  results 
with  varying  Upper  and  Lower  Bounds  for  Substructure  damping 
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C0 

^3 subs 

ub 

^3 subs 

lb 

^>sys 

ub 

sys 

lb 

^subs  max 

^3 subs 

min 

^ sys 

max 

jsys 

min 

Iterations 

Function 

Call 

0.2 

0.05 

0.01 

0.1 

0.01 

0.0500 

0.0500 

0.0498 

0.0353 

1 

42 

0.2 

0.05 

0.01 

0.04 

0.01 

0.0500 

0.0401 

0.0400 

0.0352 

2 

63 

0.2 

0.05 

0.01 

0.03 

0.01 

0.0434 

0.0301 

0.0300 

0.0298 

2 

63 

0.2 

0.05 

0.01 

0.03 

0.02 

0.0434 

0.0301 

0.0300 

0.0298 

2 

63 

0.2 

0.05 

0.01 

0.03 

0.03 

0.0434 

0.0301 

0.0300 

0.0291* 

2 

63 

0.2 

0.05 

0.01 

0.015 

0.01 

0.0217 

0.0150 

0.0150 

0.0149 

2 

63 

Table  12.  Using  a  20  Element  Beam  Substructures  and  Keeping  only  20  modes  and  Initial 
value  for  the  design  variable  kept  at  0.2;  System  and  Substructure  damping  results 
with  varying  Upper  and  Lower  Bounds  for  System  damping 


r 

jsubs 

ub 

r 

^siibs 

lb 

^ sys 

ub 

^ sys 

lb 

C  ,  max 

jsubs 

r 

^ subs 

min 

^ sys 

max 

^sys 

Min 

Iteration 

s 

Function 

Call 

1 

.015 

0.01 

0.05 

0.02 

0.0208* 

0.0208* 

0.0207 

0.0147 

2 

63 

0.2 

0.05 

0.01 

0.045 

0.04 

0.0507* 

0.0093* 

0.0456* 

0.0093* 

2 

63 

0.2 

0.05 

0.045 

0.045 

0.04 

0.0507* 

0.0443* 

0.0456* 

0.0358* 

2 

63 

0.2 

0.05 

0.01 

0.045 

0.038 

0.0500 

0.0451 

0.0450 

0.0353* 

2 

63 

0.2 

0.05 

0.045 

0.045 

0.036 

0.0500 

0.0451 

0.0450 

0.0353* 

2 

63 

0.2 

0.05 

0.045 

0.045 

0.035 

0.0500 

0.0451 

0.0450 

0.0353 

2 

63 

0.2 

0.06 

0.045 

0.045 

0.040 

0.0600 

0.0451 

0.0450 

0.0419 

2 

63 

0.2 

0.1 

0.2 

0.02 

0.05 

0.0724 

0.0501 

0.05 

0.0497 

2 

63 

Table  13.  Using  a  20  Element  Beam  Substructures  and  Keeping  only  20  modes;  System  and 
Substructure  damping  results  with  varying  Upper  and  Lower  Bounds  for  System 
and  Substructure  damping  and  Start  Point 


Size  of 
Matrix 

Time 

In  Seconds 
Trial  1 

Time  in 
Seconds 
0(n3)  Trial  1 

Time 

In  Seconds 
Trial  2 

Time  in 
Seconds 
0(n3)  Trial  2 

Time 

In  Seconds 
Trial  3 

Time  in 
Seconds 
0(n3)  Trial  3 

N 

0.000433 

0.000244 

0.000409 

0.000244 

0.000439 

0.000244 

2N 

0.002612 

0.001953 

0.002563 

0.001953 

0.002677 

0.001953 

4N 

0.017596 

0.015625 

0.017392 

0.015625 

0.017575 

0.015625 

8N 

0.127353 

0.125000 

0.126400 

0.125000 

0.128845 

0.125000 

16N 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

Table  14.  Table  showing  the  time  it  takes  to  calculate  the  eigenvectors  and  eigenvalues  of 

varying  matrix  sizes. 


Using  the  “tic,  toe,  rand,  and  eig”  commands  in  MATLAB,  we  compared  the  time 
it  takes  to  solve  for  eigvenvalues  and  eigenvectors  of  varying  size  symmetric  matrices. 
The  “tic  toe”  command  in  MATLAB,  measures  the  amount  of  time  that  it  takes 
MATLAB  to  complete  one  or  more  operations,  and  displays  the  time  in  seconds.  The 
“rand”  command  on  the  other  hand,  returns  a  pseudorandom,  scalar  value  drawn  from  a 
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uniform  distribution  on  the  unit  interval,  for  this  thesis  varying  matrix  of  sizes  N,  2N,  4N, 
8N  and  16N  comprised  of  these  random  numbers  are  compared.  The  “eig”  command 
solves  for  the  eigenvalues  and  eigenvector  of  the  given  square  matrices.  Lastly,  the  table 
also  make  use  of  the  generalized  equation  0(n  )  which  solves  for  the  amount  of  time  is 
takes  to  solve  varying  size  of  symmetric  matrices,  where  n  is  the  dimension  of  the  matrix. 
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V.  CONCLUSION 


Since  it  is  not  often  that  we  conduct  vibrations  test  on  system  structure,  the 
method  of  conducting  test  on  systems  sub  structural  parts  then  synthesizing  the  results  is 
an  efficient  way  of  predicting  or  detennining  system  modal  properties.  In  this  thesis  we 
will  exploit  the  concept  of  substructure  coupling  which  makes  use  of  both  linear 
homogenous  equations  of  constraints  and  generated  substructure  modal  parameters  such 
as  natural  frequencies,  modal  damping  and  mode  shapes  to  determine  system’s  dynamic 
response.  In  order  to  achieve  such  task,  we  made  use  of  the  zero-eigenvalue  theorem,  of 
which,  made  use  of  a  set  of  given  constraints  to  generate  a  transformation  matrix  which 
was  utilized  to  convert  modal  block  matrix,  either  of  and/or 2 Cru,  to  solve  for  system’s 
modal  parameters.  The  results  of  mentioned  coupling  are  tabulated  in  tables  one  through 
eight  of  the  results  section. 

Looking  at  tables  one  through  eight,  we  can  see  that  by  applying  more  modal 
parameters  from  each  substructure  in  solving  for  the  system  dynamic  properties,  the 
solution  converges  to  the  actual  system  solution,  which  was  generated  to  check  the  results 
from  coupling.  However,  it  is  not  often  desired  to  use  all  if  not  most  modes  in 
calculating  for  system’s  response  due  to  that  fact  that  computational  cost  becomes  a 
bigger  concern  when  solving  higher  order  matrix  problems.  Therefore,  if  we  are  able  to 
compute  the  desired  output,  in  this  case,  system’s  damping  ratios,  by  making  use  of  only 
a  few  modes  extracted  from  each  substructure,  we  are  be  able  to  eliminate  the  economic 
concerns  for  this  calculation.  Furthermore,  these  tables  also  gave  system  bounds  given 
substructure  modal  properties,  which  were  used  as  the  initial  building  block  in 
formulating  an  optimization  routine  in  solving  for  optimum  values  for  substructure 
damping  which  would  satisfy  prescribed  bounds  for  system  damping. 

Using  a  built  in  command  in  MATLAB,  “fmincon”,  and  by  slightly  modifying  the 
formulated  MATLAB  code  used  for  structural  coupling,  we  were  able  to  achieve  the  goal 
of  being  able  to  produce  values  for  substructure  damping  which  satisfies  bounds  imposed 
on  system  damping.  By  looking  at  tables  nine  through  thirteen,  which  are  results  taken 
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from  coupling  a  beam  structure  with  20  elements,  by  varying  either  start  points  or  initial 
values  for  the  design  variable  (substructure  damping)  or  the  upper  and  lower  bounds  for 
substructure  and  system  damping  ratios,  we  were  able  to  solve  for  an  optimum  value  for 
the  design  variable;  therefore  we  can  say  that  the  methodology  worked.  The  first  part  of 
the  optimization  analysis  dealt  with  varying  only  the  start  point  while  holding  the  values 
of  the  bounds  for  both  substructure  and  system  damping.  After  conducting  this  analysis 
we  may  say  that  the  start  point  did  not  have  any  impact  on  the  optimizer’s  solution 
finding  since  it  found  the  same  solution  for  all  of  the  attempts  that  was  made,  this  is  due 
to  the  fact  that  we  have  a  linear  objective  function  of  which  had  no  local  extremes  within 
the  objective  function.  Secondly,  by  looking  closely  on  tables  twelve  and  thirteen,  we 
can  see  that  the  optimizer  was  not  able  to  produce  values  that  met  the  prescribed  bounds. 
In  table  twelve  row  five  column  nine  when  we  designated  a  system  without  bounds  in 
this  case  both  lower  and  upper  bounds  were  valued  at  0.03,  the  optimizer  did  not  find  a 
feasible  solution  for  the  problem;  furthermore,  looking  at  table  thirteen  in  the  first  five 
rows  we  observe  that  values  of  either  substructure’s  or  system’s  maximum  or  minimum 
modal  damping  solutions  does  not  meet  the  bounds  set  in  the  optimization  program. 
Upon  careful  inspection  of  these  problems,  we  conclude  that  the  reason  for  not  having  a 
feasible  solution  is  because  the  bounds  set  for  these  trials  are  too  stringent  for  the 
optimizer  to  solve,  the  reason  being,  as  we  slowly  space  out  the  bounds  for  the  system 
damping  the  better  the  solution  becomes  as  shown  in  rows  six  and  seven,  row  six  showed 
how  expanding  the  bounds  for  the  system  enable  us  to  process  feasible  results  while  row 
seven  was  gave  more  separation  between  the  upper  bounds  of  system  and  substructure 
damping.  The  formulation  of  such  method  which  determines  practical  spacing  for  system 
damping  may  be  further  investigated  and  can  be  a  follow  on  topic  for  this  thesis. 

Lastly,  for  this  thesis  we  also  investigated  plate  structures  that  were  modeled  in 
PATRAN  and  NASTRAN.  During  this  process  we  were  able  to  come  up  with  modal 
parameters  for  each  substructure  after  modeling,  in  fact  we  were  able  to  convert  the  files 
generated  in  NASTRAN  to  a  useable  database  file  that  was  easily  extracted  for  use  in 
MATLAB;  however  upon  coupling  the  two  substructures  we  were  not  able  to  match 
system’s  modal  properties  which  lead  to  the  conclusion  that  there  might  be  a  mistake 
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hidden  within  the  MATLAB  code  that  was  previously  generated.  As  a  result  we  will 
continue  to  do  the  work  for  the  plate  structure;  however,  the  results  will  not  be  shown  for 
this  thesis. 
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VI.  RECOMMENDATION 


Further  investigate  why  it  is  that  system’s  maximum  natural  frequency  and 
maximum  damping  ratio  does  not  exceed  that  of  either  substructure’s  maximum 
frequency  and  damping. 

Formulate  a  method  which  determines  practical  spacing  for  system  and/or 
substructure  modal  damping  bounds  which  will  give  practical  results  that  may  be  used  in 
real  world  scenario. 

Apply  above  concept  to  more  than  two  substructures.  Create  a  multi  jointed 
structure  and  divide  them  into  substructure  then  make  use  of  the  transformation  matrix 
and  optimizer  to  produce  modal  damping  for  substructures  that  will  then  satisfy  the  given 
constraints. 
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APPENDIX 


NOTE:  If  the  main  program  calls  for  subprograms,  the  subprogram  used  for 
substructure  A  is  identical  substructure  B  with  the  exception  everything  that  refers  to  A 
must  be  replaced  with  B;  as  a  result  only  subprogram  for  substructure  A  are  shown. 


A-l  INITIAL  WORK  COUPLING  1  DOF  MASS  SPRING  SYSTEM,  WHICH 
MAKES  USE  OF  A  SYSTEM  MATRIX  THAT  WAS  CONVETED  FROM 
2ND  ORDER  DIFFERENTIAL  TO  FIRST  ORDER  DIFFERENTIAL 


%T  h  i  s  will  calculate  a  simple  2  DOF  mass  spring  system 
%S  o  I  u  t  i  on  will  be  system  natural  freq  and  damping 

%T  h  e  Solution  will  be  compared  to  the  results  gotten  from  Mathemati  ca 

%T h e  equations  of  motions  were  transformed  from  a  2nd  order  differential  to 

%a  first  order  differential 

%Gi  ven 

Ml  =  .01; 

M2  =  .02; 

K1  =  100; 

K2  =  150; 

zet  a  S u b A  =  ,02; 

zet  a'SubB  =  .01; 

%C  a  f  c  u  I  ate  Natural  Frequencies 
Wn  SubA  =  sqr  t  (  K 1  /  Ml) 

Wn'SubB  =  s  q  r  t  (  K  2  /  M2 ) 

[  Pfii  A,  LamA]  =  ei  g(  Kl,  Ml) 

Wn  A  =  s  q  r  t  (  L  a  mA ) 

[  Phi  B,  LamB]  =  ei  g(  K2,  M2) 

Wn  B  =  s  q  r  t  ( L  a  mB ) 

%C  a  I  c  u  I  ate  Cl  and  C2  using  Substructure's  zeta  and  Wn 
Cl  =  2  *  z  e  t  a  S  u  b  A  *  Wn  S  u  b  A*  Ml ; 

C2  =  2  *  z  e  t  a'SubB*Wn'SubB*M2; 

A  =  [  0,  0,  ‘ 1,  0;  07  0,  0,  1;  .  .  . 

Ml\ (  1/  150*(  Ml\ Kl)  )  -  Ml\ Kl,  -  Ml\  (  1/  1 5  0*  {  M2\  K2)  )  ,  Ml  \  (  1/  1 5 0  *  (  Ml\  C 1 )  )  -  Ml\  C 1 , 
Ml\  (  1/  1 5 0 *  (  M2 \  C2 )  )  ;  .  .  . 

-  M2\ (  1/ 150*  (  Ml\ Kl)  )  ,  M2  \  (  1/  1 5  0  *  (  M2  \  K2)  )  -  M2  \  K  2 ,  -  M2  \  {  1/ 150*{  Ml\ Cl)  )  , 

M2\  (  1/  150*  (  M2\  C 2 )  )  ■  M2\  C2] 

Lam  =  ei  g(  A)  ; 
s  i  g ma  =  real  ( Lam) ; 

Wd  =  i  ma  g  ( L  a  mj ; 
for  i  =1:1 

WD  {  i  ,  1 )  =  Wd  (  i  *2  -  1,  1)  ; 

SI  G MA ( i  ,  1)  =  si  g ma (  i  *2-  1,  1)  ; 
end 

for  i  =1:1 

Zet  a  (  i  ,  1 )  =  sqrt  ( (  (  (  WD(  i  ,  1)  A2/ SI  GMA(  i  ,  1)  A2)  +1)  A-  1)  ) ; 

WN ( i  ,  1 )  =  -  SI  G MA (  i  ,  1)  /  Zet  a(  i  ,  1) ; 
end 
WN 

Zet  a 

%Compare  with  Mathemati  ca 
Mathemati  ca  Fn  =  91.2871 

Mathemati  ca^Zeta  =  .0  1  3  6  2  7  5 
%C  heck  for  Mat  he  mat  i  ca 
C3  =  C1+C2; 

M3  =  M1+M2; 

K3  =  K 1  +  K  2 1 

Check  Fn  ='  sqrt  (  K3/  M3) 

Check'Zet  a  =  C 3 /  ( 2*  Chec  k_  F n*  M3 ) 
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A-2  COUPLING  2  DOF  MASS  SPRING  SYSTEM,  WHICH  MAKES  USE  OF  A 
SYSTEM  MATRIX  THAT  WAS  CONVETED  FROM  2nd  ORDER 
DIFFERENTIAL  TO  FIRST  ORDER  DIFFERENTIAL 


%G  i  v  e  n  Values 
nel  =2; 
nnel  =2; 
ndof  =1 ; 

n n o d e  ={  nnel  - 1)  * n e I  +1; 
sdof  =nnode*ndof ; 
conek  =  1; 
num  subs  =  2; 
ml  =  [  .  0 1 ,  .  0  2  ]  ; 
kl  =  [  1  0  0,1  2  5  ]; 
zet  asl  =  [ .  01,  .  02] ; 

m2  =  [  .  0  3 ,  ,  0  4  ]  ; 
k2  =  [  1  5  0,1  7  5  ]; 
zet  a s 2  =  [ .  01,  .  02] ; 

Kl=zeros(sdof, sdof); 

K2=zeros(sdof, sdof); 

%B u i I  d i  n g  of  K  and  M 
index  =  [0  ,  1] ; 

“/©Substructure  A 
for  i  =  1  :  length(kl); 
i  ndex  =  i  ndex  +  1; 

ka  =  [  kl(  1,  i )  ,  -  k 1 ( 1,  i  )  ;  -  k  1  (  1,  i )  ,  k 1 (  1,  i  )  ] ; 

Kl( i  ndex, i ndex)  =  Kl(i  ndex, i  ndex)  +  ka; 
end 

%A  p  p  I  y  boundary  conditions  for  subs  A 
Keep  =  2 : s  dof ; 

Ka  =  K 1  (  Keep,  Keep)  ; 

Ma  =  d  i  a  g  (  ml ) ; 

%C a  I  c  u I  a t  i  on  of  Substructure  11  A 11  s  mo d e  shapes,  natural  f  r  e q ,  and 
%p  h  y  s  i  cal  C  ma  t  r  i  x 
[  P n i  A ,  La  mA ]  =  e i  g (  K a ,  Ma )  ; 

Wa  =  s  q  r  t  (  L  a  mA )  ; 

Check_MassNor  mA  =  Phi  A1  *Ma*Phi  A; 

%S  ubstructure  A 1  s  C  ma  t  r  i  x 
for  i  =  1:  I  ength(zetasl) 

Z_Wn  A (  i  ,  i  )  =  2  *  Wa  (  i  ,  i  )  *  z  et  a  s  1  (  1 ,  i  )  ; 

end 

Ca  =  Ma*  Phi  A*  Z  Wn  A*  P  h  i  A1  *Ma; 

index  =  [0  ,  I] ; 

“/©Substructure  B 
for  i  =  1  :  I  engt  h(  k  2 )  ; 
i  ndex  =  i  ndex  +  1; 

kb  =  [  k 2 (  1,  i  )  ,  -  k 2 (  1,  i  )  ;  -  k 2 (  1,  i  )  ,  k 2 (  1,  i  )  ] ; 

K2( i  ndex, i ndex)  =  K2(i  ndex, i  ndex)  +  kb; 
end 

%A  p  p  I  y  boundary  conditions  for  subs  B 
Keep  =  1 : s  dof - 1; 

Kb  =  K 2  (  Keep,  Keep)  ; 

Mb  =  d  i  a  g  (  m2 ) ; 

%C a  I  c  u I  a t  i  on  of  Substructure  "  B "  s  mo d e  shapes,  natural  f  r  e q ,  and 
%p  h  y  s  i  cal  C  ma  t  r  i  x 
[  Phi  B,  LamB]  =  ei  g(  Kb,  Mb)  ; 

Wb  =  s  q  r  t  (  LamB)  ; 

Check  MassNor  mB  =  Phi  B1  *  Mb*  P  h  i  B; 

%S  ubstructure  B 1  s  C  ma  t  r  i  x 
for  i  =  l:length(zetas2) 

Z  Wn  B (  i  ,  i  )  =  2  *  Wb  (  i ,  i  )  *  z  e  t  a  s  2  (  1 ,  i  )  ; 

end 

Cb  =  Mb*  Phi  B*  Z  Wn  B*  P  h  i  B1  *Mb; 

Pal  =  Phi  A(  ;  ,  I)  ;  ' 

Pa  2  =  Phi  A(  ;  ,  2) ; 

Pbl  =  Phi  B{  :  ,  1)  ; 

P b 2  =  Phi  B(  :  ,  2) ; 

Mai  =  1; 

Ma  2  =  1; 

Mb  1  =  1; 

Mb  2  =  1; 

Q  =  Pal(  2,  1)  A2/ Mai  +  P  a  2 (  2,  1)  A2/ Ma  2  +  P  b 1 ( 1,  1)  A2/ Mbl  +  P  b  2 ( 1,  1)  A2/ Mb  2 ; 

%A  matrix  used  to  calculate  Wn  and  zeta  for  the  overall  structure 

%T  h  i  s  matrix  is  derived  from  the  two  substructures 

A  =  [  0,  0,  0,  0,  1,  0,  0,  0;  .  .  . 

0,  0,  0,  0,  0,  1,  0,  0;... 


%  n u mb e r  of  e I  e  me n t  s 
%  number  of  nodes  per  element 
%number  of  dofs  per  node 
%  total  number  of  nodes  in  system 
%  t  o  t  a  I  s  y  s  t  e  m  d  o  f  s 
%#  of  connection 
%#  o  f  substructure 

%ma ss  values  for  subsl 
%s  t  i  f  f  n  e  s  s  values  for  subsl 
%z  e  t  a  values  for  subsl 
%ma  ss  values  for  s  u  b  s  2 
%s  t  i  f  f  n  e  s  s  values  for  s  u  b  s  2 
%z  e  t  a  values  for  s  u  b  s  2 


ma  t  r  i  c  e  s  for  each  substructure; 
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0,  0,  0,  0,  0,  0,  1,  0;... 

0,  0,  0,  0,  0,  0,  0,  1;... 

Mal\ ( (  Q\  ( P a  1  ( 2,  1)  A2) )  *  Wa  (  1,  1)  *2)  -  Wa(l,l)A2, 

Mal\  ( (  Q\  ( P  a  1  ( 2,  1)  *  Pa2(  2,  1) ) )  *  Wa  ( 2, 2)  A2) ,  .  .  . 

-  Ma 1 \  (  (  Q\ (  P  a  1 (  2,  1) *  P  b  1  (  1,  1)  ) ) *  Wb  (  1,  1)  A2)  ,  ■ 

Mal\  (  (  Q\  (  P a  1  (  2,  1 ) *  P b 2 (  1 ,  1)  )  )  *  Wb (  2,  2)  A2)  ,  .  .  . 

Ma  1\  (  (  Q\  (  P  a  1  (  2 ,  1 )  A2 ) ) *  Z  Wn  A(  1 ,  1 )  )  -  Z  Wn  A(l,l), 

Ma  1  \  (  (  Q\  (  P  a  1  (  2 ,  1)  *  P  a  2  (  2,1)))  *Z  Wn  A(  2 ,  2 )  ) ,  .  .  7  " 

-  Ma  1  \  ( (  Q\  ( P a  1  ( 2 ,  1)  *  Pbl  ( 1"  1)  I)  *  Z  Wn  B(  1,  1) ) ,  - 

Ma 1\  (  (  Q\  (  P  a  1  (  2 ,  1 )  *  P  b  2  (  1 ,  1 ) ) )  *  Z _  Wn _  B  (  2  7  2  )  T ;  .  .  . 

Ma  2\  ( (  Q\  ( P  a  2  ( 2,  1)  *  P  a  1  ( 2,  1) ) )  *  Wa  (  1,  1)  A2) ,  Ma  2\  ( (  Q\  ( P a  2  ( 2,  1)  A2) )  *Wa(  2,  2)  A2)  - 
Wa  (  2,  2)  A2,  .  .  . 

-  Ma  2\  (  (  Q\  (  P a 2 (  2,  1)  *Pbl(  1,  1)  ) )  *  Wb (  1,  1)  A2)  ,  ■ 

Ma  2\  (  (  Q\  (  P  a  2  (  2,  1)  *  P  b  2  (  1,  1)  )  )  *Wb(  2,  2)  A2)  ,  .  .  . 

Ma  2\  (  (  Q\  (  P  a  2  (  2,  1)  *  Pa  1(  2,  1)  )  )  *Z  Wn  A(  1 ,  1 )  )  ,  Ma  2\ (  (  Q\  (  P  a  2  (  2 ,  1 )  A2 )  )  *  Z  Wn  A(  2 ,  2 )  ) 

■  Z  Wn  A(  2,  2)  ,  .  .  . 

“  Ma  2\  (  (  Q\  (  P  a  2 (  2,  1) *  P  b 1 (  1,  1) )  )  *Z  Wn  B(  1,  1) )  ,  - 

Ma  2  \  (  (  Q\  (  P  a  2  (  2 ,  1 )  *  P  b  2  (  1 ,  1 )  )  )  *  Z _  Wn  _  B  (  2  7  2  )  7;  .  .  . 

-  Mb  1  \  (  (  Q\  (  P  b  1  (  1,  1)  *  Pa  1 (  2,  1) )  )  *  Wa  (  1,  1)  A2)  ,  ■ 

Mb  1  \  ( (  Q\  (  P  b  1  (  1,  1 )  *  P  a  2  (  2 ,  1 ) )  )  *  Wa  (  2 ,  2 )  A2 )  ,  .  .  . 

Mb 1\  (  (  Q\  (  P  b 1  (  1,  1)  A2)  )  *  Wb  (  1,  1)  A2)  -  Wb  (  1 ,  1 )  A2 , 

Mb 1\  (  (  Q\  (  P  b  1  (  1,  1)  *  P  b  2  (  1,  1)  )  )  *  Wb (  2,  2)  A2)  ,  .  .  . 

-  Mb  1\  ( (  Q\  ( P b  1  (  1,  1)  *  Pa  1{  2,  1) ) )  *  Z  Wn  A(  1,  1) ) ,  - 

Mb  1\  (  (  Q\  (  P  b  1  {  1,  1)  *  P  a  2  (  2,  1)  )  )  *2  Wn  A(  2  7  2  )  T _ 

Mb  1\  (  (  Q\  (  P  b  1  (  1 ,  1 )  A2 ) )  *  Z  Wn  B(  1,  1)  )  -  Z  Wn  B(  1,  1) , 

Mb  1\  ( (  Q\  (  P  b  1  (  1,  1)  *  P b 2  (  1,  1)  )  )*Z_Wii_B(  2,  2)  ) ;  .  .  7  ' 

-  Mb2\ (  (  Q\  (  P b 2 (  1,  1)  *  Pa  1(  2,  1)  ) )  *  Wa  (  1,  1)  A2)  ,  ■ 

Mb  2\  (  (  Q\  (  P  b  2 (  1,  1)  *  Pa2(  2,  1)  ) )  *  Wa  (  2,  2)  A2)  ,  .  .  . 

Mb  2\  (  (  Q\  (  P  b  2  (  1,  1)  *  P  b  1  (  1,  1)  )  )  *  Wb  (  1,  1)  A2)  ,  Mb2\((Q\(Pb2(l,l)A2))*Wb(2,2)A2)  - 

Wb (  2,  2)  A2,  .  .  . 

-  Mb 2\  (  (  Q\  (  P b 2 (  1,  1)  *  Pa  1  (  2,  1)  )  )  *Z  Wn  A(  1 ,  1 )  )  ,  - 

Mb  2\  (  (  Q\  (  P  b  2  (  1,  1)  *  P  a  2  (  2,  1)  )  )  *2  Wn  A(  2  7  2  )  7  _ 

Mb2\((Q\(Pb2(l,  1)  *Pbl(  1,  I)  )T*Z  Wn  B(  1,  1)  )  ,  Mb  2  \  (  (  Q\  (  P  b  2  (  1 ,  1 )  A2 ) )  *  Z  Wn  B(  2,  2)  ) 
-  Z  Wn  B  (  2 ,  2 )  ]  ; 

L A M  =  ei  g(  A)  ; 

%  %T  r  i  a  I 
%  keeep  =  5:8 
%  kep  =  1:4 
% 

%  G  =  A( keeep,  keeep) 

%  L  =  A( keeep,  kep) 

%  Omega  N  =  sqrt  (  L) 

%  e i  g  (  G7 
%  ei  gj  Omega  N) 

%  f  o  r  i  =1:4 
%  Zeta  = 

%  e  i  g  e  n  v a  I  u  e s  of  A  are  s  i  g  ma  ( i  )  +/  -  j  wd  ( i  ) 

%  s  i  g  ma  =  z  e  t  a  ( i  )  *  Wn  ( i  ) 

%  Wd  ( i  )  =  Wn  (  i  j  *  s  q  r  t  ( 1  -  z  e  t  a  ( i  ) ) 

%  t  wo  eqns  in  t  wo  unknowns 
s i  gma  =  real  ( LAM) ; 

Wd  =  i  ma  g  (  LAM)  ; 

for  i  =1:4 

WD (  i  ,  1 )  =  Wd(  i  *2  -  1,  1)  ; 

SI  G MA ( i  ,  1)  =  si  g ma (  i  *2-  1,  1) ; 
end 
WD; 

SI  GMA; 

for  i  =  1:  sdof 

Zet  a  (  i  ,  1 )  =  sqrt  ( (  (  (  WD(  i  ,  1)  A2/ SI  GMA(  i  ,  1)  A2)  +1)  A-  1)  ) ; 
end 

for  i  =  1:  sdof 

WN ( i  ,  1 )  =  -  SI  G MA (  i  ,  1)  /  Zet  a(  i  ,  1) ; 
end 
WN; 

Zeta; 

Tabl  e_  S  o  I  u  t  i  on  WN  Zeta  =  [  WN  Zet  a] 

%C heck  Natural  Frequencies 
%B  u  i  I  d  i  n  g  overall  C,  M,  K  matrices 
C  =  zerosf  sdof, sdof) ; 
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M  =  zerosf  sdof, sdof) ; 

K  =  zeros! sdof, sdof) ; 
index  =  [1  ,  2]; 

C( i  ndex, i  ndex)  =  C( i  ndex, i  ndex)  +  Ca; 

M(  i  ndex,  i  ndex)  =  M(  i  ndex,  i  ndex)  +  Ma ; 

K(  i  ndex,  i  ndex)  =  K(  i  ndex,  i  ndex)  +  K a ; 

i  ndex  =  i  ndex  +  1; 

C( i  ndex, i  ndex)  =  C( i  ndex, i  ndex)  +  Cb; 

M(  i  ndex,  i  ndex)  =  M(  i  ndex,  i  ndex)  +  Mb; 

K(  i  ndex,  i  ndex)  =  K(  i  ndex,  i  ndex)  +  Kb; 

[  Phi  ,  Lam]  =  e i  g (  K ,  M) ; 
for  i  =1:  I  ength(M) 

Check  Natfreq(i.l)  =  sqrt  ( Lam(  i  ,  i  ) ) ; 
end 

Check_Natfreq; 

%  C  h  e  c  k  Z  e  t  a 
C  =  Phi  1  *C*Phi  ; 
for  i  =  1 :  I  e  n  g  t  h  (  M) 

Check  Z  e  t  a ( i ,  1 )  =  C(i , i )/(2*Check  Nat f r eq( i , 1) ) ; 
end 

Check_Zet  a; 

Tabl  e_Check_Natfreq_Zeta  =  [CheckNatfreq  CheckZet  a] 
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A-3  COUPLING  OF  2  DOF  MASS  SPRING  SYSTEM  USING  THE  ZERO- 
EIGENVALUE  THEOREM 


nel  =2;  %  number  of  el  e men t  s 

n n e I  =2 ;  %  number  of  nodes  per  element 

n d o f  =  1 ;  %number  of  dofs  per  node 

n  n  o  d  e  =(  nnel  - 1)  *  nel  +1;  %  total  number  of  nodes  in  system 

sdof=nnode*ndof;  %  total  system  dofs 

c o n e k  =  1 ;  %#  of  connecti  on 

n  u  m_  s  u  b  s  =  2 ;  %#  o  f  substructure 

ml "  =  [.01, ,02];  %ma ss  values  for  subsl 

k  1  =  [  1  0  0,1  2  5  ];  %s  t  i  f  f  n  e  s  s  values  for  subsl 

z  e  t  a  s  1  =  [  .  0 1 ,  .  0 1  ] ;  %z  e  t  a  values  for  subsl 
m2  =  [.01, ,02];  %ma  ss  values  for  s  u  b  s  2 

k  2  =  [  1  2  5,1  0  0];  %s  t  i  f  f  n  e  s  s  values  for  s  u  b  s  2 

z  e  t  a  s  2  =  [  .  0 1 ,  .  0 1  ] ;  %z  e  t  a  values  for  s  u  b  s  2 

Kl=zeros(sdof, sdof); 

K2=zeros(sdof, sdof); 

%B u i  I  d i  n g  of  K  and  M  ma t  r  i  c  e s  for  each  substructure; 
index  =  [0  ,  1]; 

%S  u  b  s  t  r  u  c  t  u  r  e  A 
for  i  =  1  :  length(kl) 
i  ndex  =  i  ndex  +  1  ; 

ka  =  [  k  1  (  1,  i )  ,  -  k 1 (  1,  i  )  ;  -  k 1 (  1,  i  )  ,  k 1 (  1,  i  )  ] ; 

Kl(i  ndex, i  ndex)  =  Kl(i  ndex, i  ndex)  +  ka; 
end 

%A p p I y  boundary  conditions  for  subs  A 
Keep  =  2 : s  dof ; 

Ka  =  K 1  (  Keep,  Keep)  ; 

Ma  =  d  i  a  g  (  ml ) ; 

%Calculation  of  Substructure  "  A "  s  mo  d  e  shapes,  natural  f  r  e q ,  and 
%p  h  y  s  i  cal  C  ma  t  r  i  x 
[  Pni  A,  LamA]  =  e  i  g  (  K  a ,  Ma )  ; 

Wa  =  s  q  r  t  (  La  mA )  ; 

Check  MassNormA  =  Phi  A1  *Ma*Phi  A; 

Checks  =  Phi  A1  *  Ka  *  P  h  i  A; 

"/(Substructure  A 1  s  C  ma  t  r  i  x 
for  i  =  1;  I  ength(zetasl) 

Z  Wn  A(  i  ,  i )  =  2*Wa(i,i)*zetasl(l,i); 

end 

Ca  =  Ma  *  P  h  i  A*Z  Wn  A*Phi  A1  *Ma; 

CheckCa  =  Phi  Ar*Ca*Phi  A; 

i  ndex  =  [ 0  ,  1] ; 

"/(Substructure  B 
for  i  =  1  :  I  engt  h(  k  2 ) 
i  ndex  =  i  ndex  +  1  ; 

kb  =  [  k 2 (  1,  i  )  ,  -  k 2 (  1,  i  )  ;  -  k 2 (  1,  i  )  ,  k 2 (  1,  i  )  ] ; 

K2( i  ndex, i ndex)  =  K2(i  ndex, i  ndex)  +  kb  ; 
end 

%Apply  boundary  conditions  for  subs  B 
Keep  =  1 :  s  dof  - 1; 

Kb  =  K 2  (  Keep,  Keep)  ; 

Mb  =  d  i  a  g  (  m2 ) ; 

“/(Calculation  of  Substructure  “  B 11  s  mo  d  e  shapes,  natural  f  r  e q ,  and 
%p  h  y  s  i  cal  C  ma  t  r  i  x 
[  Phi  B,  LamB]  =  ei  g(  Kb,  Mb)  ; 

Wb  =  sqrt(  LamB)  ; 

CheckMassNor  mB  =  Ph  i  B1  *  Mb*  P  h  i  B; 

CheckKb  =  Phi  B1  *  Kb*  P h i  B; 

“/(Substructure  B 1  s  C  ma  t  r  i  x 
for  i  =  1 :  I  engt  h(  z  et  as  2) 

Z  Wn  B (  i  ,  i  )  =  2  *  Wb  (  i  ,  i  )  *  z  e  t  a  s  2  (  1 ,  i  )  ; 

end 

Cb  =  Mb*  Phi  B*  Z  Wn  B*Phi  B1  *Mb; 

CheckCb  =  Phi  Br  *  CB  *  P  h  i  B; 

%  Solution 
Pal  =  Phi  A(  :  ,  1)  ; 

Pa  2  =  Phi  A(  :  ,  2) ; 

Pbl  =  Phi  B(  :  ,  1)  ; 

P b 2  =  Phi  B(  :  ,  2) ; 

Z  Wn _  s  u b s  =  z e r os ( s dof +1 ,  s d of +1 ) ; 

W  subs  =  zeros( sdof  +1, sdof  +1) ; 

A'=[  Pali  2,  1)  P a 2 ( 2,  1)  -  P b  1  (  1 ,  1)  ■  P b 2 ( 1 ,  1)  ] ; 

A  =  A1  *  A; 

B  =  n  u  I  I  (  A) ; 
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A*  B; 

index  =  [1,2]; 

Z  Wn  subs(  i  ndex,  i  ndex)  =  Z  Wn  subs(index,  index)  ■  Z  Wn  _  A ; 
W"s  u B s  ( i  ndex,  i  ndex)  =  W_s  u B s  ( F ndex,  i  ndex)  -  LamA; 

i  ndex  =  i  ndex  +  2; 

Z  Wn  subs(  i  ndex,  i  ndex)  =  Z  Wn  subs(  i  ndex,  index)  ■  Z  _  Wn  _  B ; 
W  subs(  i  ndex,  i  ndex)  =Wsubs(Fndex,  index)  -  LamB; 

S  =  B'*Z  Wn  sub$*B; 

L  =  B1  *W  s  u  b  S  *  B ; 

E  =  [ zeros( 3 ,  3 )  ,  e  y  e ( 3 , 3 )  ; 

L,  G]  ; 

LAM  =  e i  g (  E) 

%  e  i  g e n v a  i  u e s  of  A  are  s  i  g  ma  ( i  )  +/  -  j  wd ( i  ) 

%  s  i  g  ma  =  z  e  t  a  ( i  )  *  Wn  ( i  ) 

%  Wd  ( i  )  =  Wn  (  i  j  *  s  q  r  t  ( 1  •  z  e  t  a  ( i  ) ) 

%  t  wo  eqns  in  t  wo  unknowns 

s  i  g  ma  =  r  e  a  I  ( LAM) ; 

Wd  =  i  ma g (  LAM)  ; 
for  i  =1:3 

WD(  i  ,  1)  =  Wd(  i  *2  -  1,  1)  ; 

SI  G MA ( i  ,  1)  =  si  gma(  i  *2- 1,  1) ; 
end 
WD; 

SI  GMA; 

WD  =  sort  (  WD)  ; 

SI  GMA  =  sor  t  (  -  SI  GMA) ; 
for  i  =  1:  sdof 

Zet  a (  i  ,  1 )  =  sqrt  ( (  (  (  WD(  i  ,  1)  A2/ SI  GMA(  i  ,  1)  A2)  +1)  A-  1)  ) ; 
end 

for  i  =  1:  sdof 

WN (  i  , 1 )  =  SI  GMA( i  ,  1)  /  Zet  a(  i  ,  1)  ; 
end 
WN  ; 

Zet  a; 

Table  Solution  WN  Zet  a  =  [  WN  Zet  a] 

%C  h  e  c  k  Nat  u  r  a  F  Frequencies 
%B  u  i  I  d  i  n  g  overall  C,  M,  K  matrices 
C  =  zeros( sdof , sdof)  ; 

M  =  zerosf  sdof, sdof) ; 

K  =  zeros( sdof, sdof) ; 
index  =  [1  ,  2]; 


C(i 

ndex, i  ndex) 

=  C ( i  ndex 

i  ndex) 

+ 

Ca 

M(  i 

ndex, i  ndex) 

=  M(  i  ndex 

i  ndex) 

+ 

Ma 

K(i 

ndex, i  ndex 

=  K  ( i  ndex 

i  ndex) 

+ 

Ka 

i  ndex  =  i  ndex 

+  1 ; 

C(i 

ndex, i  ndex) 

=  C ( i  ndex 

i  ndex) 

+ 

Cb 

M(  i 

ndex, i  ndex) 

=  M(  i  ndex 

i  ndex) 

+ 

Mb 

K(i 

ndex, i  ndex) 

=  K ( i  ndex 

i  ndex) 

+ 

Kb 

[  Phi  ,  Lam]  =  e i  g (  K ,  M) ; 
for  i  =1:  I  ength(M) 

Check  Natfreq(i.l)  =  s  q  r  t  ( L  a  m(  i  ,  i  ) ) ; 
end 

Check_Nat  f  r  eq; 

%  C  h  e  c  k  Z  e  t  a 
C  =  Phi  1  * C *  P h i  ; 
for  i  =  1 :  I  e  n  g  t  h  (  M) 

Check  Zeta(i.l)  =  C(  i  ,  i  )  /  ( 2*  Check  Nat  f  r  eq(  i  ,  1) ) ; 
end' 

Check_Zet  a; 

Table  Check  Natfreq  Zeta  =  [Check  Natfreq  CheckZet  a] 

%%%%%%%Why  does  eig[G)/2Wn  =  zeta  for  overall  structure %%%% 

[  ei  gvecL,  ei  gval  L]  =  e i  g (  L ) ; 

D  i  a  g  _  K  =  ei  gvecL1  *  L  *  e  i  gvecL; 

Omega  N  =  s  o r  t  ( s  q r  t  ( d i  a g ( -  e  i  g v a  I  L) ) ) ; 

[  ei  gvecG,  ei  gval  G]  =  ei  g(  G) ; 

ANewG  =  ei  gvecL1  *G*ei  gvecL; 

Ne wG  =  sort  ( -  di  ag(  ANewG) ) ; 

%  ei  gNewG  =  ei  g(  NewG) 

%  NewGa  =  ei  gvecL'  *NewG*ei  gvecL 
%  G 

%  ei  gval  NewG  =  sort  ( -  ei  g(  NewG) ) 

%  e  i  gval  G  =  sort(-diag(G)j 
%  e  i  g  v  e  c  L ; 

for  i  =  1:  I  engt  h(  Omega  N) ; 

ei  gG  over  twoOmega  N  ( i  ,  1 )  =  NewG(  i  ,  1)  /  ( 2*Omega  N  ( i  ,  1 ) ) ; 
end 

tabl  e_  Ome  g  a  _  N_  a  n  d  _  e  i  g  G_  o  v  e  r  _  t  wo  Ome  g  a  _  N  =  [  0  me  g  a  _  N  ei  gG_over_twoOmega_N] 
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A-4  COUPLING  OF  2  ELEMENT  BEAM  STRUCTURE  USING  THE  ZERO- 
EIGENVALUE  THEOREM,  WHICH  GAVE  THE  BASIC  FOUNDATION 
FOR  THE  OPTIMIZATION  CODE 


multiply  #  elements  by  2 


system 


%#  of  e I  e me n t  s 
%#  n  o  d  e  s  /  e  I  e  me  n  t 
%#  dof/ node 
%#  total  nodes  of 
%t  o  t  a  I  s  y  s  t  e  m  d  o  f 
%L  e  n  g  t  h  of  beam  in  Inches 
%X -sectional  Area  2X2  inches 
%E  I  a  s  t  i  c  Modul  us  A I  umi  num 
"/Moment  of  Inertia  about  the  z-axis 
%E  I  e  me  n  t  Length 

%De  nsity  2702  Kg/  mA3 


%B  u  i  I  d  an  index  matrix  to  be  used  in 
ma  t  r  i  c  e  s 

R  h  o  A )  ;  %B  u  i  I  d  E I  e  me  n  t  ma  t  r  i  x 
k ,  mm,  m,  i  n  d  e x )  ;  %A s  s  e  mb  I  y  of  M 


f  o  r 
and 


and  K 


Keep  =  1:4;  %t  o  keep  all  DOF 
%l  n  p  u  t  for  Substructure  A 
%Substructure  A 
e  =2 ; 
n  p  e  =2 ; 
dof  p n  =2 ; 

t  n  =  ( e+1)  *(  npe-  1) ; 
sysdof  =  t  n *  n p e ; 

Lenght  =  10*12; 

A  =  6*6; 

E  =  1 0  e  6 ; 

I  =  1/  12*2*2A3; 

I  =  Lenght  /  e; 

RhoA  =  .0  0  0  2  8  4  9  7; 

%D  a  mp  i  n  g  Ratios 
zet  asA  =  [  .  0  0  5,  .  0  1,  .  0  1  5,  .  0  2  ]  ; 
k k  =z e r  o s ( sysdof , sysdof ) 
mm=z  er os(  sysdof ,  sysdof) 
for  i  e  I  =1 :  e 

i  n  d  e  x  =  i  ndexi  ng  s  u  b  A  ( i  el  ,  npe,  dofpn); 
building  the  overall  mass  and  stiffness 
[  k ,  m]  =b u i  I  d  MK  el  ement_subA(  E,  I  ,  I  ,  A 
( kk,  mm]  =assembf  e  MK  matrix  s  u  b  A  ( k  k 
ma  t  r  i  c  e  s 
end 

[  K  a ,  Ma  ]  =a  p  p  I  y  BC  subAfkk,  mm,  sysdof);  “/(.Applying  boundary  condition  for  the  beam 
structure 

[  Phi  A ,  La  mA  ]  =  e  i  g  (  K  a ,  Ma )  ;  “/.Calculate  natural  freq  and  mode  shapes 
Wn  A  =  s  q  r  t  (  L  a  mA )  ; 

Check  MassNorm  =  Phi  A1  *Ma*Phi  A;  “/Check  make  sure  this  is  equal  to  [I] 

Diag  R  =  Phi  A1  *Ka*Phi  A;  “/This  should  be  diagonalized  Ka  matrix 
%S  ubstructure  A 1  s  C  ma  t  r  i  x 

[  Z_  Wn  A ,  C  a  ]  =b  u  i  I  d  Ca(  zet  asA,  WnA,  Phi  A,  Ma) ;  "/.Build  physical  damping  matrix  Ca  and 
modal  damping  matrix  2*zetaA*WnA 

Check_Ca  =  Pni  A1  *Ca*Phi  A;  “/This  should  be  diagonalized  Ca  matrix 
%lnput  for  Substructure  B 
“/Substructure  B 
e  =2 ; 
n  p  e  =2 ; 
dof  p n  =2 ; 

t  n  =  ( e+1)  *(  npe-  1) ; 
sysdof  =  t  n  *  n  p  e ; 

Lenght  =  10*12; 

A  =  6*6; 

E  =  1 0  e  6 ; 

I  =  1/12*2*2 A3;  "/Moment  of  Inertia  about  the  z-axis 

I  =  Lenght  /  e; 

RhoB  =  .0  0  0  2  8  4  9  7; 

“ZD  a  mp  i  n  g  Ratios 
zet  asB  =[  .  0  0  5,  .  0  1,  .  0  1  5,  .  0  2]  ; 

kk=zeros( sysdof , sysdof ) ; 
mm=z  er  os(  sysdof ,  sysdof ) ; 
for  i  e  I  =1 :  e 

i  n  d  e  x  =  i  ndexi  ng  s  u  b  B  ( i  el  ,  npe,  dofpn); 
building  the  overall  mass  and  stiffness 
[  k ,  m]  =  b u i  I  d  MK  e I  e me nt _ s u b B(  E,  I  ,  I  ,  A, 

[  k  k ,  mm]  =a  s  s  e  mbl  e  MK  matrix  s  u  b  B  (  k  k ,  k 
end 

[  Kb,  Mb]  =appl  y  BC  subBfkk,  mm,  sysdof);  “/Applying  boundary  condition  for  the  beam 
structure 

[PhiB.LamB]  =  e  i  g  (  Kb,  Mb )  ;  “/Calculate  natural  freq  and  mode  shapes 
Wn  B  =  s  q  r  t  (  L  a  mB )  ; 

Check  MassNorm  =  Phi  B1  *Mb*Phi  B;  “/Check  make  sure  this  is 
Diag  R  =  PhiB'*Kb*PhiB;  “/This  should  be  diagonalized  Ka 
“/Substructure  B 1  s  C  ma  t  r  i  x 

[Z  Wn  B ,  C  b  ]  =b  u  i  I  d  Cb(  zet  asB,  WnB,  Phi  B,  Mb)  ;  "/Build  physical 
mo  dal  damping  matrix  2*zetaB*WnB 

Check  Cb  =  Pni  B1  *Cb*Phi  B;  “/This  should  be  diagonalized  Ca  matr 


%#  of  e I  e me n t  s 
%#  n  o  d  e  s  /  e  I  e  me  n  t 
%#  d  o  f  /  n  o  d  e 

%#  total  nodes  of  system 
"/total  s  y  s  t  e  m  d  o  f 
“/Length  of  beam  in  Inches 
%X -sectional  Area  2X2  inches 
"/Elastic  Modul  us  A I  umi  num 
“ZMo ment  of  Inertia  about  the 
%E  I  e  me  n  t  Length 

"ZDensi  t  y  2  7  0  2  Kg/  mA3 


"/Build  an  index  matrix  to  be  used  in 
ma  t  r  i  c  e  s 

RhoB);  “/Build  Element  matrix  for  M  and  K 
,  mm,  m,  i  n  d  e x )  ;  %A s  s  e  mb  I  y  of  M  and  K  ma  t  r  i  c  e  s 


equal  to  [  I  ] 
ma  t  r  i  x 

damping  matrix  Cb  and 


Sol  ut  i  on 
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[L,G,Table  Solution  WN  Z  e  t  a  ]  =  CalcWN  Zet  a(  Keep,  Phi  A,  ,  .  . 

Phi  B,  sysdof  ,  Z  Wn  A,  Z_Wn_B,  La  mA ,  L  a  mB  I ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%C heck  Natural  Frequencies %%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Building  overall  C,  M,  K  ma trices  and  solving  for  Check  Zeta  and 
%  Chec  k  Nat  f  r  eq 

[  K ,  M,  C ,  Table  Check_Natfreq_Zeta]  =bui  I  d  overall  KMC(Ca,  Ma,  Ka,  Cb,  Mb,  Kb,  sysdof,  npe,  dof 

P  n ,  e ) ; 

%Why  does  e  i  g  (  G )  /  2  Wn  =  zeta  for  overall  structure 
[ei  gvecL,  ei  gval  L]  =  eig(L); 

D i  a g _ K  =  ei  gvecL1  *L*ei  gvecL; 

Omega_N  =  sqrt  ( di  a g ( -  ei  gval  L) ) ; 

NewG  =  ei  gvecL1  *G*ei  gvecL; 

N e wG_ s  o  I  n  =  -  d  i  a  g (  NewG) ; 
for  I  =  1:  I  engt  h(  Omega  N) ; 

NewG  over  twoOmega  N  ( i  ,  1 )  =  NewG  sol  n(i  ,  l)/(2*Omega_N(i  ,  1)); 
end 

tabl  e  Omega  N  and  NewG  over  twoOmega  N  =  [Omega  N  NewG  over  twoOmega  N] 
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A-4-1  Subprogram  “indexingsubA” 


function  [  i  ndex]  =i  ndexi  ng  subA(  i  el  ,  npe,  dof  pn) 
d o f  _ e  =  npe’dof  pn; 
initial  =  ( i  el  -  1)  *(  npe- 1)  * d o f  pn; 
for  i  =1 :  d of  e 

i  n  d  e  x  (  i  )  =F n  i  t  i  a  I  +i  ; 
end 


A-4-2  Subprogram  “buildMKelementsubA” 

functi  on  [k,  m]=bui  I  d_MK_el  ement_subA(E,  I  ,  I  ,  A,  R  h  o ) 
q=E*  I  /  (  I  A3)  ; 


[  12 

6*1 

12 

6*1  . 

6*1 

4*  1 

*2  -  6*1 

2*1  A2; 

-  12 

-  6*1 

12 

-6*1;. 

6*1 

2*1 

*2  -  6*1 

4*1  ^2 ] 

=  Rho* 

A*  1  /  4  2  0 ; 

[  156 

22*1 

54 

-  13*1 

22*1 

4*1-2  13*1  -3*IA 

2;  .  .  . 

54 

13*1 

156 

-  22*1 

-  13*1  -  3*1  A2  -22*1  4*1  A2]  ; 

A  -  4  -  3  Subprogram  "assembl  e_MK  matrix  s  u  b  A " 
function  [  kk,  mm]  =a  s  s  e  mb  I  e  MKmatrix  s  u  b  A  ( k  k ,  k ,  mm,  m,  i  ndex) 
Q  =  length(index); 
for  i  =1 :  Q 
ii=index(i); 
for  j  =1 ;  Q 
jj  =i  ndex(j  ); 

kk( i i  ,  j  j  )  =kk( i i  ,  j  j  )  +k( i , j  ) ; 
end 
end 

for  i  =1:  Q 
i  i  =i  ndex(  i  ) ; 
for  j  =1:  Q 
j  i  =i  n  d  e  x ( j  )  ; 

mm(  i  i  ,  j  j  )  =mm(  i  i ,  j  j  )  +m(  i  ,  j  ) ; 
end 
end 


A-4-4  Subprogram  “applyBCsubA” 


function  [  Ka,  Ma]  =appl  y  BC  s  u  b  A  (  kk,  mm,  sysdof ) 

%Applying  Boundary  Conditions  nodes  1  and  2  are  both  zero  therefore 
%o  n  I  y  columns  and  rows  1  and  2  are  zeroes 
Keep  =  3:  sysdof ; 

Ka  =  k k ( Keep,  Keep) ; 

Ma  =  mm(  Keep,  Keep) ; 


A-4-5  Subprogram  “apply  BC  subA” 


function  [Z  Wn  A,  Ca]  =bui  I  d  _  C  a  ( zet  a  s  A ,  WnA,  Phi  A,  Ma) 
for  i  =  1:  Fengt  h(  zet  asA) 

Z_  Wn  _  A  (  i  ,  i  )  =  2*WnA(i,i)*zetasA(l,i); 

end 

Ca  =  Ma*  Phi  A*  Z_  Wn _  A*  Phi  A1  *  Ma ; 


A-4-6  Subprogram  “CalcWNZeta” 


function  [  L ,  G ,  Table  Solution  WN  Zeta]  = 

C  a  I  c  _  WN  Zeta(  Keep,  PFii  A,  Phi  B,  sysdof,  Z  Wn  A,Z  Wn  B ,  L  a  mA ,  L  a  mB ) 
%  Unco mme nt  the  following  only  if  you  want  to'compare 
%  more  than  1  number  of  modes  to  be  retained 
%  WN  =  zeros(2*(e-l)+l  ength(Ma),  length)  MODE)); 

%  Zeta  =  zeros(2*(e  -  l)+l  ength(Ma),  length)  MODE)); 

%  for  MODE 

%Wa nt  to  use  only  the  first  5  modes 
%  %  Keep  =  1:  k; 
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Z_Wn_subs  =  zeros(  2*1  ength(  Keep) ,  2*1  ength(  Keep) ) ;  %  Block  Diagonal  Matrix  of 
2*Zeta*Wn  from  Substructure  A  and  B 

W  subs  =  zeros(2*l  ength(Keep),  2*1  ength(Keep));  %B  lock  Diagonal  Matrix  of  Wn  from 
Substructure  A  and  B 

%  This  portion  will  create  an  A  matrix  (Interconnecting  Condition  Between 
%  the  Components),  the  constraint 

P h i  A ;  %Ei  gen  Vector  of  Substructure  A,  Mode  shapes 

Pal  =  Phi  A( sysdof  ■  3,:); 

P  a  2  =  Phi  A(  sysdof  ■  2,  :  )  ; 

Phi  B;  %Ei  gen  Vector  of  Su  bs  t  r  u  c  t  u  r  e  B,  Mode  shapes 

Pbl  =  Phi  B(  1,  :  ) ; 

P b 2  =  Phi  B(  2,  :  )  ; 

%  A  =  zeros(  2,  2*1  ength(  Keep)  )  ; 

A  =  [  Pal(  :  ,  Keep) , -  Pbl( : ,  Keep) ; 

P a 2 ( :  ,  Keep) , - Pb2( : ,  Keep) ] ; 

D  =  A1  *  A; 

N  =  nu  I  I  (  A) ;  %  Orthogonal  Complement  of  Dmatrix 

I  =N'  *  N;  %  Make  sure  N1  *N  =  [  I  ] 

z  e  r  o  =A*  N ;  %  Ma  ke  s  u r  e  A*  N  =  0 

%T  h  i  s  will  build  Z_Wn  subs  and  W_subs  Block  Diagonals,  inputting  values  from  each 
substructure 

i  ndex  =  [  1:  I  engt  h( Keep) ] ; 

Z  Wn  subsfindex, index)  =  Z  Wn  subs(  i  ndex,  index)  ■  Z  Wn_A(  Keep,  Keep)  ; 

W"s  uEs  ( i  ndex,  index)  =  W  suBsffndex,  index)  -  La  mA  (Keep,  Keep); 

Fndex  =  index  +  length(Keep); 

Z  Wn  subsfindex, index)  =  Z  Wn  subs(  i  ndex,  index)  ■  Z  Wn_B(  Keep,  Keep)  ; 

W"s  uEs  ( i  ndex,  index)  =  W  suBs(Tndex,index)  -  La  mB  (Keep,  Keep); 

%T  h  i  s  step  will  convert  my  diagonal  block  matrix  into  my  overall  structures 
%Mo  d  a  I  Damping  matrix  ( 2  Z  e  t  a  Wn )  and  Modal  matrix  of  eigenvalues  and  eigenvectors 
G  =  N 1  * Z_ Wn  s  u  bs  *  N;  %overall  structure 

L  =  N 1  *  W  s  u  E  s  *  N ;  %o  v  e  r  a  I  I  structure 

%  E  is  a  matrix  of  first  order  differential  equations 
E  =  [  z er os ( 2*  I  engt h ( Keep) - 2 , 2* I  engt h( Keep) - 2) ,  eye(  2*1  engt  h(  Keep)  - 
2,  2*1  engt  h (  Keep) -  2)  ; 

L,  G] ; 

LAM  =  e i g ( E ) ; 

%  e i  g e n v a  i  u e s  of  E  are  s  i  g  ma  ( i  )  +/  -  j  wd ( i  ) 

%  s  i  g  ma  =  z  e  t  a  ( i  )  *  Wn  ( i  ) 

%  Wd  (  i  )  =  Wn  (  i  j  *  s  q  r  t  ( 1  -  z  e  t  a  ( i  ) ) 

%  t  wo  eqns  and  two  unknowns 
%  Portion  Provided  by  Prof  Gordis 
si  gma  =  real  (  LAM)  ; 

Wd  =  i  ma g (  LAM)  ; 

for  i  =  1 : 2* I  engt  h( Kee  p)  -  2 
WD(  i  ,  1)  =  Wd  (  i  *2  -  1,  1)  ; 

SI  G MA ( i  ,  1)  =  si  g ma (  i  *2-  1,  1)  ; 
end 

%  WD  =  s  o  r  t  (  WD ) ;  %u  n  d  a  mp  e  d  natural  frequency 

%  SI  GMA  =  sort  (  -  SI  GMA)  ; 

%  replace  1  in  the  following  with  k-2  if  comparing 
%  mu  I  t  i  p  I  e  mo  d  e  s  Z  e  t  a  ( i  ,  1 )  and  WN  ( i  ,  1 ) 
for  i  =  1:  2*1  engt  h(  Keep)  -  2 

Zeta(i.l)  =  sqrt  ( ( ( (  WD(  i  ,  1)  *2/ SI  GMA(  i  ,  1)  A2)  +1)  A- 1) ) ; 
end 

for  i  =  1:  2*1  ength(  Keep)  -  2 

WN ( i ,  1 )  =  ■  SI  G MA ( i  ,  1 ) / Zet  a  ( i  ,  1) ; 
end 

%  Table  of  natural  frequency  and  Zeta  SOLUTION 
Tabl  e_Sol  uti  on  WNZeta  =  [  WN  Zeta] 

A  -  4  -  7  "Subprogram  Tbui  Id  overal  I  _KMC" 
f  unct i  on 

[K,M,C, Table  Check  NatfreqZeta]  =bui  I  d  overall  KMC(Ca,  Ma,  Ka,  Cb,  Mb,  Kb,  sysdof,  npe,  dof 
p  n ,  e ) ; 

C  =  zeros!  4*e-  2, 4*e-  2) ; 

M  =  zerosi  4*e-  2, 4*e-  2) ; 

K  =  zeros(  4*e-  2,  4*e-  2)  ; 
dof  e  =  e* 2; 
i  n i  t  f  a  I  =0; 
for  i  =1:  dof  e 

i  n d  e  x (  i  )  =i  ni  t  i  a  I  +i  ; 
end 

C( i  ndex, i  ndex)  =  Cj i  ndex, i  ndex)  +  Ca; 

M(  i  ndex,  i  ndex)  =  M(  i  ndex,  i  ndex)  +  Ma; 

K(  i  ndex,  i  ndex)  =  K(  i  ndex,  i  ndex)  +  Ka; 

i  ndex  =  index  +  ( 2  *  e  -  2 ) ; 

C( i  ndex, i  ndex)  =  C( i  ndex, i  ndex)  +  Cb; 

Mi  i  ndex,  i  ndex)  =  M(  i  ndex,  i  ndex)  +  Mb; 
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K(i  ndex,  i  ndex)  =  K(  i  ndex,  i  ndex)  +  Kb; 

[  Phi  ,  Lam]  =  e i  g (  K ,  M) ; 
for  i  =1:  I  engtn(M) 

Check  Natfreqfi ,  1)  =  sqrt  ( Lam(  i  ,  i ) ) ; 
end 

Check_Natfreq; 

%  C  h  e  c  k  Z  e  t  a 
C  =  Phi  1  *C*Phi  ; 
for  i  =  1 :  I  e  n  g  t  h  (  M) 

Check  Zetaji.l)  =  C( i , i ) / ( 2* Chec k  Nat f r eq( i , 1) ) ; 
end 

Check  Zet  a; 

%T a b f e  Natural  frequencies  and  Zetas  to  be  compared  to  the  solution 
Tabl  e_Check  NatfreqZeta  =  [CheckNatfreq  CheckZet  a] 


A-4-7  Subprogram  “System” 


%S e c o n d  Check  for  the  System 
e  =4 ; 
n  p  e  =2 ; 
dof  p n  =2 ; 

t  n  =  ( e  +1 ) *( npe-  1) ; 
sysdof  =  t  n *  n p e ; 


%#  of  e I  e me n t  s 
%#  n  o  d  e  s  /  e  I  e  me  n  t 
%#  dof/ node 
%#  total  nodes  of 
%t  o  t  a  I  s  y  s  t  e  m  d  o  f 


system 


Lenght  =  20*12; 

A  =  6*6; 

E  =  1 0  e  6 ; 

I  =  1/  1 2  *  2  *  2  A3 ; 

I  =  Lenght  /  e; 

RhoA  =  .0  0  0  2  8  4  9  7; 

kk=zeros( sysdof , sysdof )  ; 
mm=zeros( sysdof , sysdof ) ; 
for  i  e  I  =1 :  e 


%L  e  n  g  t  h  of  beam  in  Inches 

%X -sectional  Area  2X2  inches 

%E  I  a  s  t  i  c  Modul  us  A I  u  mi  num 

“/Moment  of  Inertia  about  the  z-axis 

%E  I  e  me  n  t  Length 

%De n s i ty  2  7  0  2  Kg/  mA3 


index  indexing  subA(iel,npe,dofpn);  %B  u  i  I  d  an  index  matrix  to  be  used 
building  the  over  a  f  I  mass  and  stiffness  ma  trices 

[  k ,  m]  =  b  u  i  I  d  MK  el  ement_subA(  E,  I  ,  I  ,  A,  RhoA)  ;  %Build  Element  matrix  for  M 
[  k  k ,  mm]  =assembf  e  MK  matrix  s  u  b  A  (  kk,  k,  mm,  m,  i  ndex) ;  “/.Assembly  of  M  and  K 
ma  t  r  i  c  e  s 
end 

Keep  =  3; sysdof ; 

Ka  =  k k ( Keep,  Keep) ; 

Ma  =  mm(  Keep,  Keep)  ; 

Keep  =  1 :  I  eng t  h(  Ka) -  2 
Kone  =  Ka( Keep,  Keep) 

Mone  =  Maj  Keep,  Keep) 


i  n 

and  K 
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A-5  DECOUPLING 


Keep  =  1:6;  %to  keep  all 
%l  n  p  u  t  for  Substructure 
%Substructure  A 
e  =3 ; 
n  p  e  =2 ; 
dot  p n  =2 ; 

t  n  =  ( e+1)  *(  npe-  1) ; 
sysdof  =  t  n *  n p e ; 

Lenght  =  10*12; 

A  =  2*2; 

E  =  1 0  e  6 ; 

I  =  1/  1 2  *  2  *  2  A3 ; 

I  =  Lenght  /  e; 

RhoA  =  .  1; 


DOF  multiply  #  elements  by  2 
A 

%#  of  e I  e me n t  s 
"/#  n  o  d  e  s  /  e  I  e  me  n  t 
"/#  dof/ node 

%#  total  nodes  of  system 

%t  o  t  a  I  s  y  s  t  e  m  d  o  f 

%L  e  n  g  t  h  of  beam  in  Inches 

%X -sectional  Area  2X2  inches 

%E  I  a  s  t  i  c  Modul  us  A I  u  mi  num 

"/Moment  of  Inertia  about  the  z-axis 

%E  I  e  me  n  t  Length 

%De  n s  i  t  y  2  7  0  2  Kg/  mA3 


kk=zeros( sysdof , sysdof )  ; 
mm=z  er  os(  sysdof ,  sysdof ) ; 
for  i  e  I  =1 :  e 


index  indexing  subA(iel,npe,dofpn);  %B  u  i  I  d  an  index  matrix  to  be  used 
building  the  overall  mass  and  stiffness  ma trices 

[  k ,  m]  =  b  u  i  I  d  MK  el  ement_subA(  E,  I  ,  I  ,  A,  RhoA);  %B  u  i  I  d  Element  matrix  for  M 
[  k  k ,  mm]  =assembf  e  MK  matrix  s  u  b  A  (  kk,  k,  mm,  m,  i  ndex) ;  "/Assembly  of  M  and  K 
ma  t  r i  c  e  s 
end 

[  Ka,  Ma]  =appl  y  BC  s  u  b  A  ( k  k ,  mm,  sysdof);  "/Applying  boundary  condition  for  the 
structure 

[  Phi  A ,  L  a  mA  ]  =  e  i  g  (  K  a ,  Ma )  ;  "/Calculate  natural  freq  and  mode  shapes 
Wn  A  =  s  q  r  t (  L  a  mA ) ; 

Check  MassNorm  =  Phi  A1  *Ma*Phi  A;  "/Check  make  sure  this  is  equal  to  [I] 

D  i  a  g  _  1?  =  Phi  A1  *Ka*Phi  A;  "/This  should  be  diagonalized  Ka  matrix 
%%%%%%%%%%%%%%%%%%!  n  p  u  t  for  Substructure  B  %%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%% S u  b  s  t  r  u  c  t  u  r  e  B  %%%%%%%%%%%%%%%%%%%%%%%%%% 


e  =3 ; 
n  p  e  =2 ; 
dof  p n  =2 ; 

t  n  =  ( e+1)  *(  npe-  1) ; 
sysdof  =  t  n *  n p e ; 

Lenght  =  10*12; 

A  =  2*2; 

E  =  1 0  e  6 ; 

I  =  1/  1 2  *  2  *  2  A3 ; 

I  =  Lenght  /  e; 

RhoB  =  .1; 

kk=zeros( sysdof , sysdof ) ; 
mm=zeros( sysdof, sysdof); 
for  i  e  I  =1 ;  e 


%#  of  e I  e me n t  s 
%#  n  o  d  e  s  /  e  I  e  me  n  t 
%#  d  o  f / n  o  d  e 

%#  total  nodes  of  system 
"/total  s  y  s  t  e  m  d  o  f 
"/Length  of  beam  in  Inches 
%X -sectional  Area  2X2  inches 
"/Elastic  Modul  us  A I  u  mi  num 
“/Moment  of  Inertia  about  the  z-axis 
%E  I  e  me  n  t  Length 
"/Density  2  7  0  2  Kg/mA3 


index  indexing  subB(iel,npe,dofpn);  "/Build  an  index  matrix  to  be  used  i 
building  the  overall  mass  and  stiffness  ma  trices 

[  k ,  m]  =  b  u  i  i  d  MK  el  ement_subB(  E,  I  ,  I  ,  A,  RhoB)  ;  "/Build  Element  matrix  for  M 
[  kk,  mm]  =assemble  MK  matrix  s  u  b  B  (  kk,  k,  mm,  m,  i  ndex) ;  "/Assembly  of  M  and  K 
end 

[  Kb,  Mb]  =appl  y  BC  s  u  b  B  ( k  k ,  mm,  sysdof);  "/Applying  boundary  condition  for  the 
structure 

[  Phi  B ,  L  a  mB  ]  =  e  i  g  (  Kb,  Mb )  ;  "/Calculate  natural  freq  and  mode  shapes 
Wn  B  =  s  q  r  t  (  L  a  mB ) ; 

Check  MassNorm  =  Phi  B1  *Mb*Phi  B;  "/Check  make  sure  this  is  equal  to  [I] 

Diag  R  =  PhiB'*Kb*PhiB;  "/This  should  be  diagonalized  Ka  matrix 
%S uBstructure  B's  C  matrix 

%  [Z  Wn  B ,  C  b  ]  =  b  u  i  I  d  CbfzetasB,  WnB,  Phi  B,  Mb);  "/Build  physical  damping  matrix 
modal _ d a mp i  n g  matri  x_2*zetaB*WnB 

%  Check  Cb  =  Phi  B'  *Cb*Phi  B;  "/This  should  be  diagonalized  Ca  matrix 
%  Uncomment  the  folio  wing  only  if  you  want  to  comp  are 
%  mo  re  than  1  number  of  modes  to  be  retained 
%  WN  =  zeros(  2*(e  -  1)  +1  ength(  Ma)  ,  length)  MODE)); 

%  Zeta  =  zeros!  2*(  e  -  1)  +1  engt  h(  Ma) ,  length)  MODE)); 

%  for  MODE 


%Wa nt  to  use  only  the  first  5  modes 
%  %  Keep  =  1;  k; 

W  subs  =  zeros!  2*1  ength(  Keep) ,  2*1  engt  h(  Keep) ) ;  "/Block  Diagonal  Matrix  of 
SuBst  r  uct  ur  e  A  and  B 

"/This  portion  will  create  an  A  matrix  (Interconnecting  Condition  Between 
"/the  Components),  the  constraint 

PhiA;  "/Eigen  Vector  of  Substructure  A,  Mode  shapes 

Pal  =  Phi  A(  sysdof  -  3,;); 

P  a  2  =  Phi  A(  sysdof  -  2,  :  )  ; 


i  n 

and  K 

beam 


n 

and  K 
matri  c  e  s 

beam 


Cb  and 


Wn  from 
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Phi  B;  %Ei  gen  Vector  of  Substructure  B,  Mode  shapes 

Pbl  =  Phi  B(  1,  :  ) ; 

P b 2  =  Phi  B(  2, :  ) ; 

A  =  [  P a  1 ( : ,  Keep) , - P b  1  ( : ,  Keep) ; 

P a 2 ( : ,  Keep) , - P b 2 ( : ,  Keep) ] ; 

D  =  A1  *A 

N  =  nu  I  I  (  D)  %  Orthogonal  Comp  I  ement  of  D  mat  r  i  x 

I  n  =N 1  *  N  %  Make  sure  N1  *  N  =  [  I  ] 

zero=D*N  %  Make  sure  A*N  =  0 

%T  h  i  s  will  build  Wsubs  Block  Diagonals,  inputting  values  from  each  substructure 
i  ndex  =  [ 1:  I  engt  h( Keep) ] ; 

Wsubs(index,index)  =  W  subs(index.index)  -  L  a  mA  (  Keep,  Keep)  ; 
i  ndex  =  index  +  I  engt  h(  Keep)  ; 

Wsubsfindex, index)  =  W  subsfindex,  index)  -  La  mB  (Keep,  Keep); 

%%This  step  will  convert  my  diagonal  block  matrix  into  my  system  structures 
L  =  N'*W_subs*N;  %system  structure 

%%%%%%%%%f  i  n  d  natural  f  r  e  q  for  s  y  s  t  e  m %%%%%%%%%%%% 

[  ei  gvecL,  ei  gval  L]  =  ei  g(  L) 

WN  =  sortfdi  ag(sqrt(-ei  gval  L))) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%  Sol  uti  on  %%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%G  e  n  e  r  a  t  e  Values  for  system  Zeta 
a  =  .  03; 
b  =  .  01; 
modes  =  10; 

Zeta  =  a  +  ( b-a) .  *rand(  modes,  1) ; 

%B  u  i  I  d  s  y  s  t  e  m  d  a  mp  i  n  g  ma  t  r  i  x 
for  i  =  1 ;  mo  d  e  s 

TWO  ZETA  WN  SYSTEM!  i  ,  i  )  =  2*  Zet  a  (  i  ,  1 )  *  WN(  i  ,  1 )  ; 
end 

TWO  ZETA  WN  SYSTEM  =  eigvecLHWO  ZETA  WN  SYSTEM*ei  gvecL1 
%F  i  n  d  substructures  Modal  damping  matrices 
Two  Zeta  Wn  A  B  =  N*TWO  ZETA  WN  SYSTE M* N1  ; 

WnA  =  di  ag(sort(di  ag(WnA),  'descend')); 

WnB  =  di  ag(sort(di  ag(WnA),  'descend')); 

%lgnoring  the  off  diagonals  assume  proportionally  damped  substructures 
Asubs  =z  e  r  o  s  (  6 ,  6 )  ;  Bsubs=zeros(6,  6) 
i  ndex  =  [ 1:  I  engt  h( Keep) ] ; 

Asubs(  i  ndex,  i  ndex)  =  Two  Zeta  Wn  A  B(  i  ndex,  i  ndex) 
save  =  [7;length(Keep)*2j 

Bsubs(  index,  index)  =  Two_Zeta  Wn_A  B(  7:  12,  7:  12) 

%  Ra  =  Phi  A'  *As ubs* Phi  A" 

%  Rb  =  Phi  B'  *  Bs  u bs  *  P h i  B 
for  i  =  1 :  I  e  n  g  t  h  (  Keep) 

z  e  t  a  s  A ( i ,  i  )  =  As ubs ( i , i  )  /  ( 2  *  Wn  A ( i , i  ) ) 

z  e  t  a  s  B  (  i  ,  i  )  =  Bsubs(i,i)/(2*  Wn  B ( i , i  ) ) 

end 
Zet  a 

Zet  asA  =  di  ag(zetasA) 

Zet  asB  =  di  a g ( zet  asB) 

%  [  Phi  A ,  L  a  mA  ]  =  e  i  g  (  Ka ,  Ma ) ;  "/(.Calculate  natural  freq  and  mode  shapes 
%  WnA  =  s  q  r  t  (  La  mA )  ; 

%  [  Phi  B ,  L  a  mB  j  =  e  i  g  (  K b ,  Mb ) ;  "/(Calculate  natural  freq  and  mode  shapes 
%  Wn  B  =  s  q  r  t  (  L  a  mB ) 
for  i  =  1:  I  engt  h(  zet  asA) 

Z  Wn  A(  i  ,  i  )  =  2  *  WnA ( i  ,  i  ) *  z  e  t  a  s  A(  i  ,  i  ) ; 

Z_  Wn _  B ( i ,  i  )  =  2*  Wn  B (  i  ,  i  )  * zet as B( i  ,  i  )  ; 

end 

%Check  if  it  is  a  I  mo  s  t  close  to  system  damping 
Z  Wn  subs  =  zeros!  2*1  ength(  Keep) ,  2*1  ength(  Keep) ) ; 
i  ndex”=  [  1;  I  engt  h(  Keep)  ] ; 

Z  Wn_subs(  index,  index)  =  Z  Wn  subs(  i  ndex,  index)  +  Z  Wn_A(  Keep,  Keep)  ; 
i  ndex  =  i  ndex  +  I engt  h( Keep) ; ~ 

Z  Wn  subs!  i  ndex, i  ndex)  =  Z  Wn  subs!  i  ndex, i  ndex)  +  Z  Wn  B(  Keep,  Keep) ; 

G  ='N'  *Z  Wn  subs*N; 
sol  =  di  a g (  GJ 
TWO  ZETA  WN  SYSTEM; 
check  =  df  a  g[  TWO_  ZETA_  WN_  SYSTEM) 

The  subprograms  used  for  this  decoupling  program  is  the  same  as  the  one  used  for  the  beam 

coupling. 
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A-6  FMINCON  BEAM 


global  modes  e  Uppersys  Lowersys 


%R  h  o 


Lenght  =  10*12; 
A  =  6*6; 

E  =  1 0  e  6 ; 

I  =  1/  12*2*2A3; 
Lenght / e; 

.0  0  0  2  8  4  9  7; 


%L  e  n  g  t  h  of  beam  in  Inches 
%X  -  s  e  c  t  i  o  n  a  I  Area  6X6  inches 
%EI  ast  i  c  Mo  du  I  us  A I  umi  num 
%Moment  of  Inertia  about  the  z-axis 
%E  I  e  me  n  t  Length 
%Dens i  t  y 


%  di  sp( 1  ZetaA  &  ZetaB  (damping  ratios  for  each  substructure)’) 

%  d  i  s  p  (  1  F  i  n  d  bounds  for  ZetaA  &  ZetaB  that  would  meet  the  prescribed  constraint  for 
the  System  damping’) 

%  e  =  input)’ Enter  number  of  elements  for  each  Substructure:’) 
e  =  20 


%  dispf’lnputing  desired  damping  ratios  for  each  substructure1) 

%  init  damping  =  i  nput(’  Enter  Initial  value  for  substructure  damping  typically 
bet  ween  .  02  and  .  2: 1  ) 
init  damping  =.2 

%  dfsp(’  Upper  and  Lower  bound  for  each  substructure’) 

%  Uppersubs  =  i  nput(  1  Enter  Upper  bound  for  substructure  damping:1) 

%  Lowersubs  =  i  nput(’  Enter  Lower  bound  for  substructure  damping:1) 

Uppersubs  =  .060 
Lower  s  ubs  =  .04  5 

%  dispf’Prescribed  bounds  for  system  damping1) 

%  Uppersys  =  i  nput(’  Enter  Upper  bound  for  sys  damping:1) 

%  Lowersys  =  i  nput(’  Enter  Lower  bound  for  sys  d a mp i  n g :  ’  ) 

Uppersys  =  ,045 
Lowersys  =  .040 


%  select  =  input)  ’To  use  ail  modes  enter  1.  If  not  enter  2:’) 

% 

%  i  f  select  ==  1; 

%  modes  =  e  *  2 


%  else 


%  displ’lnputing  desired  modes  kept’) 

%  modes  =  inputf’lnput  #  of  modes  to  be  Kept:1) 
%  end 

modes  =  2  0 

zetas  init  =  zeros)  modes, 1); 
for  i  =  1 :  mo  d  e  s 

zetas  i  n  i  t  ( i  ,  1 )  =  init_damping; 
end 

zet  asA  =  zetas  init; 


A=[  ] ; 
b=[  1 ; 

A e  q  =[  ]  ; 
b  e  q  =[  ] ; 
fori  = 

I  b  ( i 
u  b  (  i 
end 

[  f ,  f  v  a  I 
zet  asA 
zet  asB  = 


1 :  modes 

,  1)  =  Lower  s  ubs ; 

,  1)  =  Uppersubs; 

,  exi  tfl  ag,  output] 
=  f  ’  ; 
f  ’  ; 


fmi  neon)  @o  b  j  fun,  zetasA,  A,  b,  Aeq,  beq,  I  b,  ub,  @c  ons ) 


%lnput  for  Substructure  A 
%Substructure  A 
Keep  =  1:  modes; 
n  p  e  =2 ; 
dof  p n  =2 ; 

t  n  =  ( e+1)  *(  npe-  1) ; 
sysdof  =  t  n*npe; 

Lenght  =  10*12; 

A  =  6*6; 

E  =  1 0  e  6 ; 

I  =  1/  12  *  2*  2  A3 ; 

I  =  Lenght  /  e; 

RhoA  =  .0  0  0  2  8  4  9  7; 
k k  =z  e r  o s  ( sysdof ,  sysdof ) ; 
mm=zeros( sysdof, sysdof); 
for  i  e  I  =1 :  e 


%#  n  o  d  e  s  /  e !  e  me  n  t 
%#  d  o  f / n  o  d  e 

%#  total  nodes  of  system 
%t  o  t  a  I  s  y  s  t  e  m  d  o  f 
%L  e  n  g  t  h  of  beam  in  Inches 
%X- sect i  onal  Area  2X2  i  nches 
%EI  ast  i  c  Modul  us  A I  umi  num 
%Moment  of  Inertia  about  the  z-axis 
%E  I  e  me  n  t  Length 
%De  n  s i  t  y 


index  indexing  subA(iel,npe,dofpn);  %B  u  i  I  d  an  index  matrix  to  be  used  in 

building  the  over  a  F I  mass  and  stiffness  ma  trices 

[  k ,  m]  =  b  u  i  I  d  MK  el  ement  _s  u  b  A  (  E,  I  ,  I  ,  A,  RhoA) ;  %Bui  I  d  Element  matrix  for  M  and 


K 
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[  k  k ,  mm]  =assembl  eMKmatri  x  s  u  b  A  ( kk,  k,  mm,  m,  i  ndex) ;  %A  s  s  e  mb  I  y  of  M  and  K 
ma  t  r i  c  e  s 
end 

[  Ka,  Ma]  =appl  y_BC  s  u  b  A  ( k  k ,  mm,  s  y  s  d  o  f ) ;  “/(.Applying  boundary  condition  for  the  beam 
structure 

[  Phi  A ,  L  a  mA  ]  =  e  i  g  (  K  a ,  Ma )  ;  “/(.Calculate  natural  freq  and  mode  shapes 
Wn  A  =  s  q  r  t (  L  a  mA ) ; 

Check_MassNorm  =  Phi  A1  *Ma*Phi  A;  %C  heck  make  sure  this  is  equal  to  [I] 

D  i  a  g  _  K  =  Phi  A1  *Ka*Phi  A;  %T  h  i  s  should  be  diagonalized  Ka  matrix 


“/(.Substructure  A 1  s  C  ma  t  r  i  x 
if  modes  ==  2*e 

[Z  Wn_A,  Ca]  =bui  I  d_Ca(zetasA,  WnA,  Phi  A,  Ma)  ;  %B  u  i  I  d  physical 
modal  damping  matrix  2*zetaA*WnA 

Check  Ca  =  Phi  A1  *Ca*Phi  A;  %T  h  i  s  should  be  diagonalized  Ca 
else' 


[Z  Wn  A  ]  =  c  o  n  s  Z_Wn  A(  zet  asA,  WnA) ;  %Bui  I  d  physical 
damping  matrix  2*  zet  a  A*  WnA 


end 

%lnput  for  Substructure  B 
n  p  e  =2 ; 
dof  p n  =2 ; 

t  n  =  ( e+1)  *(  npe-  1) ; 
sysdof  =  t  n *  n p e ; 


%#  n  o  d  e  s  /  e  I  e  me  n  t 
%#  d  o  f / n  o  d  e 

%#  total  nodes  of  system 
%total  system  dof 


d  a  mp  i  n  g 


dampi  ng  mat  r  i  x  Ca  and 
ma  t  r  i  x 

matrix  Ca  and  modal 


Lenght  =  10*12; 

A  =  6*6; 

E  =  1 0  e  6 ; 

I  =  1/  12*2*2A3; 

I  =  Lenght  /  e; 

RhoB  =  .0  0  0  2  8  4  9  7; 
kk=zeros( sysdof ,  sysdof )  ; 
mm=z  er  os(sysdof,sysdof); 
for  i  el  =1:  e 


%Length  of  beam  in  Inches 
%X -sectional  Area  2X2  inches 
"/(Elastic  Mo  d  u  I  us  A I  umi  num 
"/(Moment  of  Inertia  about  the  z-axis 
%E  I  e  me  n  t  Length 
%Densi  ty 


index  indexing  subB(iel,npe,dofpn);  %B  u  i  I  d  an  index  matrix  to  be  used  in 
building  the  overall  mass  and  stiffness  ma  trices 

[  k,  m]  =bui  I  d_MK  element  subB(  E,  I  ,  I  ,  A,  RhoB) ;  %B  u  i  I  d  Element  matrix  for  M  and  K 
[  k  k ,  mm]  =assembf  e  MK  matrix  s  u  b  B  ( k  k ,  k ,  mm,  m,  i  ndex)  ;  “/.Assembly  of  M  and  K  matrices 
end 

[  Kb,  Mb]  =appl  y_BC  s  u  b  B  ( k  k ,  mm,  sysdof);  “/(Applying  boundary  condition  for  the  beam 
structure 

[  Phi  B ,  L  a  mB  ]  =  e  i  g  (  Kb,  Mb )  ;  "/(Calculate  natural  freq  and  mode  shapes 
Wn  B  =  s  q  r  t  (  L  a  mB ) ; 

Check  Ma  s  s  No  r  m  =  Phi  B1  *Mb*Phi  B;  “/Check  make  sure  this  is  equal  to  [I] 

Diag  R  =  PhiB'*Kb*PhiB;  “/This  should  be  diagonalized  Ka  matrix 
%SuBstructure  B's  C  matrix 
if  modes  ==  2*e 

[Z  Wn  B ,  C  b  ]  =  b  u  i  I  d  CbfzetasB,  WnB,  Phi  B,  Mb);  “/Build  physical  damping  matrix  Cb  and 
modal  damping  matrix  2*zetaB*WnB 

Check  Cb  =  Phi  B1  *Cb*Phi  B;  “/This  should  be  diagonalized  Ca  matrix 

else 


[Z  Wn_  B]  =c  ons_Z_Wn_  B(  z  et  as  B,  Wn  B) ;  “/Build  physical  damping  matrix  Cb  and  modal 
damping  matrix  2*zetaB*WnB 
end 


%%%%“/%%%%%  Solution  %%“/%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Z_Wn_subs  =  z  er  os  ( 2*  I  engt  h(  Keep) ,  2*  I  engt  h(  Keep) ) ;  %  Block  Diagonal  Matrix  of 
2*Zeta*Wn  from  Substructure  A  and  B 

W  subs  =  zeros(2*l  ength(Keep),  2*1  ength(Keep));  “/Block  Diagonal  Matrix  of  Wn  from 
Substructure  A  and  B 

“/This  portion  will  create  an  A  matrix  (Interconnecting  Condition  Between 
“/the  Components),  the  constraint 

PhiA;  “/Eigen  Vector  of  Substructure  A,  Mode  shapes 

Pal  =  Phi  A( sysdof  ■  3,;); 

P  a  2  =  Phi  A(  sysdof  ■  2,  :  )  ; 

Phi  B;  “/Eigen  Vector  of  Su  bs  t  r  u  c  t  u  r  e  B,  Mode  shapes 

Pbl  =  Phi  B(  1,  :  ) ; 

P b 2  =  Phi  B(  2,  :  )  ; 

%  A  =  zeros(  2,  2*1  ength(  Keep) ) ; 

A  =  [  Pal(  ;  ,  Keep) , - Pbl( ; ,  Keep) ; 

P a 2 (  :  ,  Keep) , - Pb2( : ,  Keep)  ]  ; 

D  =  A1  *  A; 

N  =  nu  I  I  (  A) ;  “/Orthogonal  Complement  of  Dmatrix 

I  =N'  *  N;  %  Make  sure  N1  *N  =  [  I  ] 

zero  =A*  N ;  %  Make  sure  A*N  =  0 

“/This  will  build  Z_Wn  subs  and  W_subs  Block  Diagonals,  inputting  values  from  each 
substructure 

i  ndex  =  [  1:  I  engt  h( Keep) ] ; 
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Z_Wn  s  ubs  ( i  ndex,  i  ndex)  =  Z  Wn  subs(  i  ndex,  i  ndex)  ■  Z_Wn_A(  Keep,  Keep)  ; 

Wsubsfindex,  index)  =  W  subsfTndex.index)  -  L  a  mA  (  Keep,  Keep)  ; 

Tndex  =  index  +  I  engt  h[ Keep) ; 

Z_Wn  s  ubs  ( i  ndex,  i  ndex)  =  Z  Wn  s  ubs  ( i  ndex,  i  ndex)  ■  Z_Wn_B(  Keep,  Keep)  ; 

Wsubsfindex,  index)  =Wsubs(Tndex,  index)  -  L  a  mB  (  Keep,  Keep)  ; 

%T  h  i  s  step  will  convert  my  diagonal  block  matrix  into  my  overall  structures 
%Mo  d  a  I  Damping  matrix  ( 2  Z  e  t  a  Wn )  and  Modal  matrix  of  eigenvalues  and  eigenvectors 
G  =  N 1  *  Z  Wn  s  u  b  s  *  N ;  %o  v  e  r  a  I  I  structure 
L  =  N '  *  W~  s  u  E  s  *  N ;  %o  v  e  r  a  I  I  structure 

"/.solving  using  a  L  and  G%%%%%  METHOD  3  IS  FASTER  THAN  METHOD  2 
[ei  gvecL,  ei  gval  L]  =  eig(L); 

Di  a g _ K  =  ei  gvecL1  *L*ei  gvecL; 

WN  =  s  q  r  t  ( d  i  a  g ( -  e  i  g  v  a  I  L ) ) ; 

NewG  =  ei  gvecL1  *G*ei  gvecL; 

N e wG_ s  o  I  n  =  -  d  i  a  g (  NewG) ; 

for  i  =  1 :  I  e  n  g  t  h  (  WN ) ; 

Zetasys(i.l)  =  NewG  sol  n(i  ,  l)/(2*  WN  ( i  ,  1 ) ) ; 
end 

Table  Solution  =  [WN,  Zetasys]; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%C heck  Natural  Frequencies %%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if  modes  ==  2*e 

%  Building  overall  C,  M,  K  matrices  and  solving  for  Check  Zeta  and 
%  Check  Natfreq 

[  K,  M,  C,  Table  Check]  =bui  I  d  overall  KMC(  Ca,  Ma,  Ka,  Cb,  Mb,  Kb,  sysdof ,  npe,  dof  pn,  e)  ; 

else 

NOCHECK  =  1 
end 

zetas  initmax  =  max(zetas_init) 
zetas~initmin  =  minjzetas  init) 
z  e  t  a  s  A  ma  x  =  ma  x  ( f ) 
zet  as  A"  mi  n  =  mi  n  ( f ) 

Zet  asys_max  =  max(Zetasys) 

Zetasys" mi  n  =  min(Zetasys) 


A-6-1  Subprogram  “Cons” 


function  [c,ceq]  =  cons(zetasA) 
global  modes  e  Uppersys  Lowersys 
zet  asA  =  zet  a s A 1 ; 
z  et  as  B  =  z  et  as  A; 

Keep  =  1:  modes; 


n  p  e  =2 ; 
dof  p n  =2 ; 

t  n  =  ( e+1)  *(  npe-  1) ; 
sysdof  =  t  n *  n p e ; 

Lengbt  =  10*12; 

A  =  6*6; 

E  =  1 0  e  6 ; 

I  =  1/ 1 2  *  2  *  2  A3 ; 

I  =  Lenght  /  e; 

RhoA  =  .00028497; 


nput  for  Substructure 

%#  n  o  d  e  s  /  e  I  e  me  n  t 
%#  d  o  f  /  n  o  d  e 

%#  total  nodes  of  system 
%t  o  t  a  I  s  y  s  t  e  m  d  o  f 
%L  e  n  g  t  h  of  beam  in  Inches 
%X -sectional  Area  2X2  inches 
%E  I  a  s  t  i  c  Modul  us  A I  umi  num 
"/.Moment  of  Inertia  about  the  z- 
%E  I  e  me  n  t  Length 
"/.Densi  ty 


kk=zeros( sysdof , sysdof )  ; 
mm=zeros( sysdof, sysdof); 
for  i  el  =1:  e 


index  indexing  subA(iel,npe,dofpn);  “/.Build  an  index  matrix  to  be  used  in 
building  the  over  a  f  I  mass  and  stiffness  ma  trices 

[  k,  m]  =bui  I  d_MK  el  ement  _s  u  b  A  (  E,  I  ,  I  ,  A,  RhoA) ;  %B  u  i  I  d  Element  matrix  for  M  and  K 
[  k  k ,  mm]  =assembf  e  MK  matrix  s  u  b  A  ( k  k ,  k ,  mm,  m,  i  ndex)  ;  "/.Assembly  of  M  and  K 
ma  t  r  i  c  e  s 
end 

[  Ka,  Ma]  =appl  y  BC  s  u  b  A  ( k  k ,  mm,  sysdof);  "/Applying  boundary  condition  for  the  beam 
structure 

[  Phi  A ,  L  a  mA  ]  =  e  i  g  (  K  a ,  Ma )  ;  "/.Calculate  natural  freq  and  mode  shapes 
Wn  A  =  s  q  r  t  (  L  a  mA )  ; 

Check_MassNorm  =  Phi  A1  *Ma*Phi  A;  "/Check  make  sure  this  is  equal  to  [I] 

Di  ag  K  =  P  h  i  A1  *  Ka*  Ph  i  A;  "/This  should  be  diagonalized  Ka  matrix 

“/Substructure  A 1  s  C  ma  t  r  i  x 

[Z  Wn  A  ]  =c  o  n  s  Z  Wn_A(  zet  asA,  WnA) ;  “/Build  physical  damping  matrix  Ca  and  modal 
damping  mat  r  i  x“2*zet  aA*WnA 

%%%%%%%%%%%%%%%%%%!  nput  for  Substructure  B  %%%%%%%%%%%%%%%% 
n p e  =2 ;  %#  nodes  /  el  ement 
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dof  p n  =2 ; 

t  n  =  ( e+1)  *(  npe-  1) ; 
sysdof  =  t  n * n p e ; 

Lenght  =  10*12; 

A  =  6*6; 

E  =  1 0  e  6 ; 

I  =  1/  12  *  2*  2  *3 ; 

I  =  Lenght  /  e; 

Rh  o  B  =  .0  0  0  2  8  4  9  7; 
k k  =z e r  o s ( sysdof , sysdof ) ; 
mm=z er os( sysdof ,  sysdof )  ; 
for  i  e  I  =1 :  e 


%#  d  o  f / n  o  d  e 

%#  total  nodes  of  system 
%t  o  t  a  I  system  dof 
%L  e  n  g  t  h  of  beam  in  Inches 
%X -sectional  Area  2X2  inches 
%E  I  a  s  t  i  c  Mo  d  u  I  us  A I  umi  num 
"/(.Moment  of  Inertia  about  the  z-axis 
%E  I  e  me  n  t  Length 
%De  n  s  i  ty 


index  indexing  subB(  i  el  ,  npe,  dofpn) ;  %B  u  i  I  d  an  index  matrix  to  be  used  in 
building  the  overall  mass  and  stiffness  ma trices 

[  k ,  m]  =b  u  i  I  d  MK  el  ement  _s  u  b  B  (  E,  I  ,  I  ,  A,  RhoB) ;  %B  u  i  I  d  Element  matrix  for  M  and  K 
[  k  k ,  mm]  =assembf  e_MK  matrix  s  u  b  B  (  kk,  k,  mm,  m,  i  ndex) ;  %A  s  s  e  mb  I  y  of  M  and  K  matrices 
end 

[  Kb,  Mb]  =appl  y  BC  s  u  b  B  ( k  k ,  mm,  sysdof);  "/(.Applying  boundary  condition  for  the  beam 
structure 

[  Phi  B ,  L  a  mB  ]  =  e  i  g  (  Kb,  Mb )  ;  "/.Calculate  natural  freq  and  mode  shapes 
Wn  B  =  s  q  r  t  (  L  a  mB ) ; 

Check  Mas  s  Nor  m  =  Phi  B1  *  Mb*  Phi  B;  %C  heck  make  sure  this  is  equal  to  [  I  ] 

Diag  R  =  PhiB'*Kb*PhiB;  %This  should  be  diagonalized  Ka  matrix 
%S  uBst  ruct  ure  B 1  s  C  ma  t  r  i  x 

[  Z_  Wn  B]=cons_Z  Wn_B(  zet  asB,  WnB) ;  %B  u  i  I  d  physical  damping  matrix  Cb  and  modal 
damping  matri  x~25zetaB*WnB 

%%%%%%%%%%  Solution  %%%%%%%%%%%% 

Z  Wn  _  s  u  b  s  =  zeros!  2*1  engt  h(  Keep)  ,  2*1  engt  h(  Keep) ) ;  %  Block  Diagonal  Matrix  of 
2*Zeta*Wn  from  Substructure  A  and  B 

W  subs  =  zeros(2*l  ength(Keep),  2*1  ength(Keep));  %B  lock  Diagonal  Matrix  of  Wn  from 
Substructure  A  and  B 

%  This  portion  will  create  an  A  matrix  (Interconnecting  Condition  Between 
%  the  Components),  the  constraint 

PhiA;  %Ei  gen  Vector  of  Substructure  A,  Mode  shapes 

Pal  =  Phi  A( sysdof  ■  3,;); 

P  a  2  =  Phi  A(  sysdof  ■  2,  :  ) ; 

Phi  B;  %Ei  gen  Vector  of  Su  bs  t  r  u  c  t  u  r  e  B,  Mode  shapes 

Pbl  =  Phi  B(  1,  ;  ) ; 

P b 2  =  Phi  B(  2,  :  )  ; 

%  A  =  zeros(  2,  2*1  ength(  Keep)  )  ; 

A  =  [  Pal(  ;  ,  Keep) , - Pbl( ; ,  Keep) ; 

Pa 2  (  :  ,  Keep) , - P b 2 ( : ,  Keep) ] ; 

D  =  A1  *  A; 

N  =  n  u  I  I  (  A) ;  %  Orthogonal  Complement  of  Dmatrix 

I  =N'  * N;  %  Make  sure  N1  *N  =  [  I  ] 

zero  =A*  N ;  %  Make  sure  A*N  =  0 

%This  will  build  Z_Wn  subs  and  W_subs  Block  Diagonals,  inputting  values  from  each 
substructure 

i  ndex  =  [  1:  I  engt  h( Keep) ] ; 

Z  Wn  subs(  i  ndex,  index)  =  Z  Wn  subs(  index,  index)  ■  Z  Wn_A(  Keep,  Keep)  ; 

W~suEs(  index,  index)  =WsuEs(Tndex,  index)  -  La  mA  (Keep,  Keep); 

Index  =  index  +  I  engt  hf  Keep) ; 

Z  Wn  subsjindex,  index)  =  Z  Wn  subs(  i  ndex,  index)  ■  Z  Wn_B(  Keep,  Keep)  ; 

W~suEs(  index,  index)  =  WsuEs(Tndex,  index)  -  L  a  mB  (  Keep,  Keep)  ; 

%Thi  s  step  will  convert  my  diagonal  block  matrix  into  my  overall  structures 
%Modal  Damping  matrix  ( 2  Z  e  t  a  Wn )  and  Modal  matrix  of  eigenvalues  and  eigenvectors 
G  =  N'*Z  Wn  s  u  b  s  *  N ;  "/(.overall  structure 
L  =  N'*W"suEs*N;  "/(.overall  structure 

%%%%s  o  I  v  i  n  g  using  a  L  and  G%%%%%  METHOD  3  FASTER  THAN  METHOD  2 
[  ei  gvecL,  ei  gval  L]  =  ei  g(  L) ; 

D i  a g _ K  =  ei  gvecL1  * L * e i  gvecL; 

WN  =  sqrtfdi  ag(-ei  gval  L)); 

NewG  =  ei  gvecL1  *G*ei  gvecL; 

N e wG_ s  o I  n  =  -  di  a g (  NewG) ; 
f  o  r  ~  i  =  l;length(WN); 

Zetasys(i.l)  =  NewG  sol  n(i  ,  l)/(2*  WN  ( i  ,  1 )  )  ; 
end 

lb  sys  =  Lowersys; 
ubsys  =  Uppersys; 

Zet  aq  =  Zet  as ys ; 

Zet  ar  =  Zet  asys ; 
c  =  z  er  os ( I  engt  h(  Zet  as  ys ) *  2, 1) ; 
for  i  =  1 :  I  engt  h(  Zet  as  ys )  *  2; 
if  i  <  I  ength(Zetasys)  +1; 

c(  i  ,  1)  =  Zet  aq(  i  ,  1)  -  u b _ s  y s ; 

else 


c ( i ,  1 )  =  ■ Zet a r  ( I  e ngt h (  Zet a s y s ) ,  1 )  +  I  b _ s  y s ; 
end 


end 
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c  e  q  =  [  ] ; 


A-6-2  Subprogram  “objfun” 

function  [f,fval  ,  exi  tfl  ag,  output]  =  objfun(x) 
global  mo  d  e  s 
f  =  0; 

for  i  =  1 :  mo  d  e  s 
f  =  f  ■  x  (  i  ) ; 
end 


A-6-3  Subprogram  “cons_Z_Wn_A” 

function  [Z  Wn  A  ]  =b  u  i  I  d  Ca(  zet  asA,  WnA) 
for  i  =  1:  Fengt  h(  zet  a  s  A ) 

Z _ Wn _ A (  i  ,  i  )  =  2*WnA(i  ,  i  )*zetasA(l,  i  ); 

end 

The  subprograms  used  for  this  FMINCONJBEAM  program  is  the  same  as  the  one  used  for 
the  beam  coupling. 
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A-7 


COMAPARING  THE  LOADWORK  IN  CALCULATING 
EIGENVALUES  OF  MATRIX  SIZES  N,  2N,  4N,  8N  AND  16N 


%C  ode  Wr  i  1 1  en  By  Prof  Go  r  d  i  s 

clear 

cl  c 

n  =  50; 

a  =  rand(n,n); 

a  =  (a  +  a ' ) / 2 ;  %  create  sy mme tic  a 

ti  c;  [phi  ,  I  am]  =  eig(a);  tl  =  toe; 

n  =  2  *  n; 
a  =  rand(n,n); 

a  =  ( a  +  a '  )  /  2 ;  %  c  r  e  a  t  e  s  y  mme  t  i  c  a 

ti  c;  [phi  ,  I  am]  =  eig(a);  1 2  =  toe; 

n  =  2  *  n; 
a  =  r  a  n  d  ( n ,  n ) ; 

a  =  ( a  +  a '  )  /  2 ;  %  create  sy  mme  tic  a 

ti  c;  [phi  ,  I  am]  =  eig(a);  1 3  =  toe; 

n  =  2  *  n; 
a  =  rand(n,n); 

a  =  (a  +  a ' ) / 2 ;  %  create  sy  mme  tic  a 

ti  c;  [phi  ,  I  am]  =  eig(a);  t  4  =  toe; 

n  =  2  *  n; 
a  =  rand(n,n); 

a  =  ( a  +  a '  )  /  2 ;  %  c  r  e  a  t  e  s  y  mme  t  i  c  a 

ti  c;  [phi  ,  I  am]  =  eig(a);  1 5  =  toe; 

1 1  =  [  t  1  t  2  t  3  t  4  t  5] 
tt  =  1 1 / 1 1 ( 5 ) ; 


nn  =  [nA3  (2*n)A3  (4*n)A3  ( 8*  n )  A3  ( 1 6*  n )  A3  ] 

nn  =  n  n / n  n ( 5 ) ; 

[  1 1 ;  n  n  ] 


THE 
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